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FOREWORD 


This  report  is  prepared  as  part  of  an  In-house  effort 
under  Project  2401,  Task  No.  240102,  "Design  and  Analysis 
Methods  for  Aerospace  Vehicle  Structures,"  and  Work  Unit 
24010208,  "Automated  Design  of  Advanced  Aerospace  Structures." 
The  work  was  carried  out  in  the  Design  and  Analysis  Methods 
Group  of  the  Analysis  & Optimization  Branch  (FBR),  Structural 
Mechanics  Division,  Air  Force  Flight  Dynamics  Laboratory, 
Wright- Patterson  AFB,  Ohio. 

The  time  period  of  the  effort  for  this  Technical  Report 


is  May  1978  - August  1978..  The  manuscript  was  originally  \ 

released  as  AFFDL- FBR- TM- 78-89  in  August  1978. 
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1 .  INTRODUCTION 


The  program  "ANALYZE"  was  orginally  developed  for  in-house  studies 
in  structural  analysis  and  optimization.  It  is  the  basis  for  a number  of 
programs  such  as  "OPTSTAT"^ , OPTCOMP^  and  "DANALYZ"^ . "ANALYZE" 
has  been  used  by  the  authors  in  their  consultation  work  on  a number  of 
Air  Force  projects.  It  has  also  been  used  as  a demonstration  program 
in  structural  analysis  courses  at  the  University  of  Dayton  and  the  Air 
Force  Institute  of  Technology.  This  program  was  distributed  earlier 
with  makeshift  input  and  output  instructions.  These  instructions  did 
not  include  details  of  the  theory  nor  the  internal  organization  of  the 
program.  The  purpose  of  this  report  is  to  generate  comprehensive  docu- 
mentation for  the  "ANALYZE"  program. 

The  program  is  based  on  the  displacement  method  of  finite  element 
analysis . In  such  an  analysis  the  continuum  is  replaced  by  a 
discrete  model  consisting  of  a finite  number  of  nodes  connected  by 
elements  (See  Figure  1).  This  discretization  reduces  the  original 
differential  equations  of  the  continuum  to  a set  of  algebraic  equations 
which  can  be  solved  much  more  readily  on  digital  computers. 

The  program  has  basically  four  finite  elements: 

1.  Bar  (Axial  Force  Member) 

2.  Membrane  Triangle 

3.  Membrane  Quadrilateral 

4.  Shear  Panel 

four  elements  and  their  local  coordinate  systems  are  shown  in  Figure  2. 
bar  is  a constant  strain  line  element  and  is  equivalent  to  a rod 

1 


The 

The 


element  in  the  NASTRAN^  program.  The  membrane  triangle  is  a constant 

strain  plate  element  similar  to  TRMEM  in  NASTRAN.  The  membrane 

quadrilateral  is  constructed  out  of  four  (non-overlapping)  constant 

strain  membrane  triangles  (element  2)  with  a fictitious  interior  node. 

This  interior  node  is  later  removed  by  static  condensation.  This  element 

is  similar  to  QDMEM2  in  NASTRAN.  The  shear  panel  is  also  constructed 

out  of  four  non-overlapping  triangles  with  a fictitious  interior  node. 

However,  only  the  shear  energy  is  considered  in  determining  the  stiffness 

of  this  element.  Although  the  formulation  is  somewhat  different,  this 

element  gives  comparable  results  to  the  NASTRAN  SHEAR  element  or  the 

(8) 

so  called  Garvey  shear  paner  1 . 

The  basis  for  the  derivation  of  the  shear  panel  is  empirical,  and  it 
is  primarily  intended  to  eliminate  some  of  the  difficulties  encountered 
in  using  membrane  triangles  and  quadrilaterals.  For  example,  in  beam 
problems  (rectangular  beams,  I-beam,  Box  Beams  including  multicell  wings 
and  fuselage  structures)  the  high  stress  gradients  in  the  webs  do  not 
justify  the  use  of  constant  strain  triangles  or  quadrilaterals  derived 
from  these  triangles.  In  fact,  use  of  such  elements  for  the  webs 
(spars  and  ribs  in  wings)  overestimates  the  stiffness  by  an  order  of 
magnitude.  Aerospace  engineers  have  offset  this  difficulty  to  a large 
extent  by  judicious  use  of  membrane  elements  in  conjunction  with  the 
shear  panels.  In  fact  the  early  finite  element  models  of  wings  and 
fuselages  consisted  primarily  of  bars  and  shear  panels.  However,  the 
present  practice  of  using  membrane  triangles  and  quadrilaterals  for  the 
top  and  bottom  skins,  bars  for  the  posts,  spar  and  rib  caps,  and  shear 


4 


panels  for  the  spars  and  ribs  eliminates  to  a large  extent  the  need 
for  determining  the  equivalent  thicknesses  and  cross-sectional 
areas  in  the  bars  and  shear  panels  model.  The  models  consisting  of 
these  elements  are  most  satisfactory  for  determining  the  primary  load 
paths  in  built-up  structures  such  as  wings  and  fuselages.  In  addition 
the  simplicity  of  these  elements  makes  interpretation  of  the  results 
easy  and  also  keeps  the  analysis  costs  low  because  the  stiffness 
matrices  of  these  elements  can  be  generated  in  a fraction  of  a second. 

The  detailed  formulation  and  additional  information  on  these  elements  are 
given  in  Section  3. 

In  finite  element  analysis,  a large  proportion  of  the  time  is 
spent  in  the  solution  of  the  force  displacement  relations.  The  program 
uses  standard  Gaussian  elimination  with  modifications  to  take  into  account 
the  symmetry  and  sparseness  characteristics  of  the  stiffness  matrix. 

The  details  of  the  solution  scheme  and  storage  of  the  stiffness  matrix 
are  given  in  Sections  2 and  5.  "ANALYZE"  is  an  incore  program  whose 
core  requirements  depend  on  the  problem  size,  primarily  measured  in  terms 
of  the  number  of  degrees  of  freedom  and  the  size  of  the  semi-bandwidth. 
However,  the  bandwidth  per  se  is  not  considered  in  the  program.  With 
an  available  core  of  about  lOOKg  one  can  solve  problems  of  up  to  300 
to  400  degrees  of  freedom.  With  the  full  core  of  a machine  like  the 
CDC  6600,  it  is  possible  to  solve  problems  of  up  to  1500  degrees  of  freedom 
and  a comparable  number  of  elements.  The  details  of  core  requirements 
are  discussed  in  Appendix  A. 

The  program  is  written  in  standard  ANSI  Fortran  IV. 
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2.  ANALYSIS 


In  the  finite  element  analysis  the  continuum  is  replaced  by  a 
discrete  model  consisting  of  a finite  number  of  nodes  connected  by 
elements  (members).  The  rationale  in  such  an  approximation  is  that  the 
response  between  the  nodes  (i.e.,  in  the  elements)  can  be  expressed  as 
a function  of  the  response  at  the  nodes.  The  functional  relationship 
between  the  two  responses  is  approximated  by  various  interpolation 
functions  or  shape  functions.  The  type  of  functions  depends  on  the 
complexity  of  the  problem  at  hand.  This  discretization  reduces  the 
original  differential  equations  of  the  continuum  to  a set  of  algebraic 
equations  which  can  be  solved  much  more  readily  on  digital  computers. 

The  equations  of  the  finite  element  analysis  can  be  derived 
conveniently  by  considering  the  strain  energy  of  the  deformed  system. 

For  example,  if  the  elastic  body  is  idealized  by  m finite  elements 
connecting  q nodes  (See  Figure  1),  the  strain  energy  of  the  ith  element 
can  be  written  as 

Ti'?v  2**5)  dv  Cl 

where  and  e..  are  the  stress  and  strain  vectors  and  Vi  is  the  volume  of 
the  element.  For  a linearly  elastic  body  the  relation  between  stress  and 
strain  can  be  written  as 


!i 


(2) 


* Superscript  t on  a matrix  represents  transpose 


where  E.  is  the  symmetric  matrix  of  material  elastic  constants.  For 
typical  plane  stress  problems  the  elastic  constants  matrix  is  of  dimension 
3x3.  For  an  isotropic  material  in  plane  stress  problems  the  elements 
of  E are  as  follows: 


E = 


1 y 0 

y 1 0 

0 0 J-O-y) 


(3) 


where  E and  y are  the  elastic  modulus  and  poisson's  ratio  of  the  material 
respectively.  For  an  orthotropic  material  the  elastic  constants  matrix  is 
given  by 

1 ye 

E = j yB  B 

~ 1-By‘ 

0 0 

where  E^  and  E?  are  the  longitudinal  and  transverse  moduli,  respectively,  in 
the  directions  of  the  material  property  axes.  B is  the  ratio  of  transverse 
to  longitudinal  modulus  ^/E-j).  G and  y are  the  shear  modulus  and 
poisson's  ratio  respectively. 

The  essence  of  the  finite  element  approximation  is  that  the  internal 
displacements  of  the  elements  are  expressed  as  functions  of  the  displacements 
of  the  discrete  nodes  to  which  they  are  connected.  The  local  coordinate 
systems  and  the  nodal  degrees  of  freedom  of  the  four  elements  are  shown 
in  Figure  2.  The  functional  relationship  between  the  element  internal 
displacements  and  the  discrete  nodal  displacements  is  given  by 

!?1  = !^i  (5) 


0 

§{1 -By2) 
L1 


(4) 
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where  the  matrix  represents  the  displacements  in  the  element  which  are 
functions  of  the  spatial  coordinates  (x,  y).  The  shape  function  <J>.j  is 
a rectangular  matrix,  and  its  elements  are  also  functions  of  the  spatial 
coordinates.  The  vector  vi  represents  the  nodal  displacements  in  the 
direction  of  the  element  degrees  of  freedom  in  the  local  coordinate 
system  (Figure  2).  Now  the  strain-displacement  relations  can  be 
written  as 

£i  = B Wi  (6) 

where  B is  a differential  operator.  For  a plane  stress  problem  B is 
given  by 

B = 


3 

3x 


0 


0 


d_ 

9y 


a. 

9y 


a_ 

ax 


(7) 


Substitution  of  Equations  2,  5 and  6 in  1 gives  the  expression  for  strain 
energy  in  the  following  form 


where  ki  is  the  element  (member)  stiffness  matrix  with  respect  to  the 
discrete  coordinates  v and  is  given  by 


(8) 


~i  = 


B1  E.  B 4>.  dV 


(9) 


An  alternate  but  a convenient  method  of  determining  the  elements  of  the 

(9) 

member  stiffness  matrix  is  by  invoking  the  principle  of  virtual  work'  7 
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which  gives 


(10) 


Jvi 

where  o|p^  is  the  stress  state  due  to  the  element  displacement  configuration 
in  which  vp  = 1 while  all  other  v's  are  zero.  Similarly  is  the  strain 
state  due  to  the  unit  displacement  configuration  in  the  direction  of  the  qth 
degree  of  freedom.  These  two  conditions  are  shown  in  Figure  3 for  the 
degrees  of  freedom  1 and  2 of  the  membrane  triangle.  It  should  be  noted 
that  besides  assuming  appropriate  shape  functions,  the  integration  in 
Equations  9 or  10  is  one  of  the  difficult  tasks  in  the  case  of  complex 
elements  in  finite  element  analysis.  However,  for  membrane  elements  this 
integration  does  not  present  any  difficulties  as  will  be  seen  in  the 
next  section.  For  more  complex  elements  the  usual  practice  is  to  adopt 
numerical  integration  schemes.^10,11  ^ 

From  Equation  8 and  Castigliano's  first  theorem,  the  relation  between 
the  element  nodal  forces  and  the  displacements  may  be  written  as 


!i 


= ki  !i 


(ii) 


where  s.  is  the  element  nodal  force  matrix  corresponding  to  the  displacement 
matrix  v^.  Similar  force-displacement  relations  for  the  total  structure 
can  be  derived  from  the  strain  energy  of  the  structure.  The  total  strain 
energy  T of  the  structure  can  be  written  as  the  sum  of  the  energies  of 
the  individual  components. 

m » m a 

r • l T,  = 2 i vi  ki  V1  (lz) 

i=l  1 L 1=1  -1  -1 
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(a)  First  Unit  Mode 


(b)  Second  Unit  Mode 


FIG.  3:  Examples  of  Unit  Displacement  Modes 


FIG.  4:  Quadrilateral  or  Shear  Panel 
Divided  into  Four  Triangles 


In  general,  for  most  structures,  it  is  convenient  to  define  a local 
coordinate  system  for  each  element  and  a global  coordinate  system  for  the 
total  structure.  In  such  a case  the  element  and  structure  generalized 
coordinates  can  be  related  by 

Vi  = a.  u (13) 

where  ai  is  the  compatibility  matrix.  Its  elements  can  be  determined 
by  kinematic  reasoning  alone  provided  the  structure  is  kinematically 
determinate.  The  matrix  u Is  the  generalized  displacement  vector  of 
the  structure  in  the  global  coordinate  system.  It  is  interesting  to 
note  that  Equation  13  not  only  transforms  element  displacements  from  local 
to  global  coordinates  but  also  gives  information  about  how  the  elements 
are  connected  to  the  structure.  From  Equation  13  and  the  principle  of 
virtual  work  it  is  easy  to  show  that  the  transformation  between  the 
forces  on  the  structure  and  the  element  internal  forces  is  given  by 

P = a*  s.  (14) 

where  P is  the  force  vector  on  the  structure  in  the  global  coordinate  system. 
The  transformation  given  in  Equation  14  is  sometimes  referred  to  as  a 
contragradient  transformation. 

Substitution  of  Equation  13  in  12  gives  the  expression  for  the  total 
strain  energy  in  the  form 

r ■ \ u*  K u (15) 


where  K,  the  total  stiffness  matrix  of  the  structure,  is  written  as  the 
sum  of  the  component  stiffness  matrices. 


Again  using  Castigliano's  first  theorem  the  relation  between  the  generalized 
force  matrix  P corresponding  to  the  displacement  matrix  u may  be  written  as 

= K u (17) 

In  most  structural  analysis  problems  the  stiffness  matrix  K is 
sparsely  populated.  It  is  essential  to  take  advantage  of  this  fact 
in  solving  the  load  deflection  equations  (Equation  17),  particularly  in 
the  case  of  problems  with  a large  number  of  degrees  of  freedom  where  the 
cost  of  computation  can  be  prohibitive  otherwise.  The  "ANALYZE"  program 
uses  Gaussian  elimination  with  modifications  to  take  into  account  the 
symmetry  and  sparseness  of  the  stiffness  matrix. 

Basically  Gaussian  elimination  involves  decomposition  of  the  stiffness 
matrix  by 

K = L D Ll  (18) 

where  L is  the  unit  lower  triangular  matrix  and  D is  a diagonal  matrix. 

The  advantage  of  this  decomposition  scheme  is  that  the  L matrix  retains 
some  of  the  sparseness  characteristics  of  K which  consequently  reduces  the 
numbers  of  computations.  Also  L and  D can  be  assigned  the  same  storage 
as  K. 

The  next  step  is  the  forward  substitution  by 

L Y = P (19) 

where  the  matrix  Y is  given  by 

Y = D u (20) 

In  Equation  19  the  solution  of  Y can  be  accomplished  by  simple  forward 
substitution.  Once  Y is  obtained,  u can  be  solved  by  back  substitution 
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7.  Determination  of  the  element  displacements  in  their  local 


l 

coordinate  system  (Equation  13). 

8.  Determination  of  the  stresses  in  each  element  (Equations  6,  5,  and  2). 

9.  Output  the  structure  displacements,  element  stresses  and  other 
information  such  as  element  strain  energies,  etc. 

The  next  section  consists  of  the  details  of  the  stiffness  matrix 
formulations  for  the  four  elements  in  this  program.  \ 
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where  and  x2  are  the  coordinates  of  the  two  ends  in  the  local  coordinate 
system.  Then  the  shape  function  (Equation  5)  corresponding  to  this  linear 
displacement  field  can  be  written  as 


[x1  - x2J 


(x  - 


x2),  -(x  - x1 


From  the  strain-disDlacement  relations,  the  axial  strain  in  the  element  is 


given  by 


_ 3w  „ , 
ex  3x 


From  the  principle  of  virtual  work  (Equation  10)  the  individual  elements  of 
the  member  stiffness  matrix  can  be  written  as 

k,j  - f o'0  e‘J)  dV  = « (21) 


where  A is  the  cross-sectional  area,  L is  the  length  of  the  member,  and  E is 
the  modulus  of  elasticity  of  the  material.  The  member  stiffness  matrix  is 


given  by 


k _ AE  1 -1 
£ “ L -1  1 


The  member  force  matrix  is  given  by 
s = k v 

The  stress  in  the  member  is  given  by 


°x  = E ex 


„ . $1  . "S2 

° " a — r 


>>vA, 


The  strain  energy  in  the  element  is  given  by 

vh‘?  <26> 

or 

Ti  = \ °x  £x  A L (27) 

TRIANGULAR  MEMBRANE  ELEMENT 

The  membrane  triangle  is  the  basic  plate  element  in  the  program.  It 
is  used  to  construct  the  membrane  quadrilateral  as  well  as  the  shear  panel 
with  some  modifications.  The  membrane  triangle  can  be  used  effectively  in 
all  cases  where  the  primary  loading  is  inplane  forces.  These  include  top 
and  bottom  skins  of  aircraft  wings,  flanges  of  I and  box  beams  when  they 
are  subjected  to  constant  normal  stresses  (tension  or  compression) 
only  and  skins  of  sandwich  construction.  However,  they  are  not  suitable 
for  situations  where  high  stress  gradients  exist.  For  example,  they  are 
unsuitable  for  spars  and  ribs  of  wings  and  other  lifting  surfaces,  webs 
of  I and  box  beams  and  flat  plates  where  the  primary  load  is  bending. 

If  used  in  such  cases,  they  overestimate  the  stiffness  or  generate  singularity. 
Figure  2 shows  the  triangle  elements  with  the  local  coordinate  system. 

The  generalized  coordinates  v^ , — » vg  ^present  the  inplane 

displacements  of  the  three  nodes  in  the  local  coordinate  system.  The 
displacement  field  in  the  element  is  assumed  to  be  linear.  This  gives 
constant  strain  in  the  element.  For  a linearly  elastic  material  the 
stress  in  the  element  will  also  be  constant. 


17 


The  linear  displacement  field  in  the  element  can  be  represented 
wx  = al  x + bl  y + C1 


by 


(28) 


wy  = a2  x + b2  y + c2 

where  wx  and  wy  are  the  x-y  displacements  in  the  plane  of  the  plate 
in  the  local  coordinate  system.  a1 , b-|  etc.  are  the  six  undetermined 
coefficients.  Equation  28  can  be  written  in  matrix  form  as  follows: 


w = 


x y 1 0 0 0 
0 0 0 x y 1 


(29) 


The  six  unknown  coefficients  can  be  uniquely  determined  by  the  six 
boundary  conditions  at  the  nodes. 
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(30) 


where  x-j,  y1 , , x3  and  y3  are  the  coordinates  of  the  three  nodes  of 

the  triangle  in  the  local  coordinate  system.  It  should  be  noted  that  the 


nodal  displacements  are  grouped  into  x and  y directions,  so  that  the 

nodal  coordinate  matrix  on  the  right  hand  side  partitions  into  a diagonal 

matrix.  The  inversion  of  the  partitioned  diagonal  matrix  involves  simply 

the  inversion  of  the  component  matrix.  Now  the  shape  matrix  <p  is  given  by 
,-l 

(31) 


= x 1 


where  the  matrix  x is  given  by 
[7y  1 0 0 01 


x = 


0 0 0 x y lj 


and  the  Z matrix  is  given  by 


Z = 


O' 


0 1 X 

I ~ 


(32) 


(33) 


The  coordinate  matrix  X is  given  by 


X = 


-1 


(34) 


It  is  interesting  to  note  that  each  column  of  Z represents  a unit 

displacement  mode:  i.e.  the  jth  column  of  the  inverse  represents  a 

displacement  mode  in  which  v.  = 1 while  all  other  nodal  displacements  are 

J 

zero  (See  Figure  3).  This  fact  is  used  to  advantage  in  determining  the 
elements  of  the  member  stiffness  matrix. 

From  linear  strain-displacement  relations  the  strains  can  be  written 
as 


f x.  3x 


= a 


(35) 


3wy 

c s — J-  = K 

y 3y  i 


(36) 


(37) 


T 


t 


9WX  9wy 

exy  = 9y~  + 3x  = ^1  + a2 

From  the  principle  of  virtual  work  (Equation  10)  the  elements  of  the 
member  stiffness  matrix  can  be  written  as 

,-\t 


k. . = 
ij 


dv  . 


V 


E ^ dV 


(38) 


where  and  are  the  stress  and  strain  matrices  corresponding  to 
the  unit  displacement  modes  explained  under  Equation  34.  Since  the  linear 
displacement  variation  implies  constant  strain,  the  integral  in  Equation  38 
can  be  replaced  by  the  volume  of  the  element: 

1 * .(IT  r Jj) 


kjj  ■ j|X|  tjW'ir 


(39) 


where  |X|  is  the  determinant  of  the  nodal  coordinate  matrix  which  represents 
twice  the  area  of  the  element  and  t is  the  thickness  of  the  element.  Now 
the  stiffness  matrix  of  the  element  is  given  by 

vt 


k = J l?l  1 


.0)' 

A2)1 

.(6)1 


i£0) 

...  el'!*;. 
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e(2)‘  e e(2) 

...  c(2)‘  E «, 

E e 
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.(6)T 


E E 


(2). 


.(6)X 


E e 


(6) 

(6) 

(6) 


(40) 


The  member  force  matrix  is  given  by 
s = k v 

The  stress  matrix  in  the  element  is  given  by 
a = E e 


(41) 


(42) 
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The  strain  energy  in  the  element  is  given  by 


Ti  = \ v (43) 

or 

Ti  = 1 l?l  t ^ 5 (44) 

The  next  important  step  in  the  evaluation  of  the  stress  state  in 

an  element  is  the  selection  of  a suitable  failure  criteria  because  of  the 

combined  stresses  (a  , a .and  a ) in  plate  elements.  The  modified  energy  of 

x y xy 

distortion  or  Von-Mises  criteria  is  adopted  to  determine  the  effective 
stress  in  an  element.  The  effective  stress  is  given  by 


aeff  = (ax  + ay  ' ax  ay  + 3axy)V2 


(45) 


The  margin  of  safety  is  evaluated  by  first  determining  the  effective  stress 
ratio  (ESR)  , ,, 


ESR  = 


a 2 a 2 a a a 2 

* ti>  - * <#> 


(46) 


where  XX  and  YY  are  the  tension  or  compression  allowable  in  the  x and  y 
directions,  respectively,  and  ZZ  is  the  shear  allowable.  Then  the  margin  of 
safety  (MS)  is  determined  by 

MS  - (47) 

If  the  user  does  not  provide  the  allowable  stress  values,  then  default 
values  of  60,000  psi  for  tension  and  compression  allowables  In  both 
directions  and  36,000  psi  for  the  shear  allowable  are  used. 

QUADRILATERAL  MEMBRANE  ELEMENT 

The  quadrilateral  element  is  most  frequently  used  to  represent  membrane 
skins  unless  the  corners  etc.  require  the  use  of  the  triangular  element. 


21 


I 


Figure  4 shows  the  local  coordinate  system  and  the  generalized  coordinates 


(displacements)  v1  through  vg.  The  element  is  assumed  to  be  a flat  plate, 
and  all  nodes  are  assumed  to  lie  on  a plane  connecting  the  first  three 
nodes  (1,  2,  and  3).  In  effect  the  warping  in  the  element  is  ignored. 

This  approximation  results  in  an  overestimation  of  the  stiffness  of  a 
truly  warped  quadrilateral  element.  In  most  cases  the  effect  of  the 
approximation  is  small,  and  it  can  be  further  reduced  by  reducinq  the  mesh 
size  of  the  model  in  the  regions  of  high  warping.  However,  if  the  warp 


is  too  large,  the  quadrilateral  should  be  broken  up  into  two  or  more 
triangles. 

As  mentioned  earlier,  the  stiffness  of  the  quadrilateral  element  is 

determined  by  breaking  it  into  four  component  triangles  as  shown  in 

Figure  4.  A fictitious  node  in  the  quadrilateral  is  located  by  averaging 

the  coordinates  of  the  four  nodes  as  given  by 

_ X,  + x2  + x3  ♦ x4  (48) 


The  stiffness  of  the  four  triangles  is  then  computed  by  Equation  40  in  the 
local  coordinate  system  shown  in  Figure  2c.  Addition  of  the  four  stiffness 
matrices  gives  a 10  x 10  stiffness  matrix  with  two  degrees  of  freedom 
included  for  the  fifth  node.  This  fictitious  node  is  later  removed  by 
static  condensation  before  adding  to  the  total  structure.  The  procedure 
for  static  condensation  is  outlined  next. 


The  force  displacement  relations  of  the  5 node  quadrilateral  are 


written  as 


Bq  = Bq  Tq 


(50) 


where  the  subscript  refers  to  the  quadrilateral  element  with  5 nodes. 
Equation  50,  partitioned  to  isolate  the  degrees  of  freedom  of  the  fifth 
node,  can  be  written  as 


5i 

h,_i_j_h2_n_ 

Ci 

Ell 

hi,  ij  hi,  n 

hi 

(51) 


Equation  51  can  be  written  as  two  separate  equations 

h = h,I  Cl  + h.II  Til 


(52) 


hi  = ~II ,1  Cl  + l5lI.II  ClI 


(53) 


Since  the  fifth  node  does  not  actually  exist  in  the  original  model,  no 
external  forces  can  be  applied  to  this  node.  This  condition  gives 


ClI  = “~II ,11  hi, I Cl 

Substitution  of  Equation  54  in  52  gives 

h = (h,i  " h,n  hr, ii  hi.i*  Ci 


(54) 


(55) 


From  Equation  55  the  stiffness  matrix  of  the  original  quadrilateral  can 
be  written  as 

B “ hi  " h,ll  hi, II  hi, I (56) 

The  stiffness  as  obtained  by  Equation  56  is  added  to  the  total  structure 
after  appropriate  coordinate  transformations  to  the  global  coordinate  system. 
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When  the  structure  displacements  are  determined,  the  fifth  node  displacements 
can  be  determined  by  Equation  54.  Now  the  stresses  in  each  triangle  can  be 
determined  as  before.  The  effective  stress  ratio  is  determined  for  each 
triangle  separately  (Equation  46),  and  then  a weighted  average  is  used 
in  computing  the  effective  stress  ratio  and  the  margin  of  safety.  This 
weighted  average  is  computed  by 

(ESRh  A,  + (ESR)2  A?  + (ESR),  A,  + (ESR),  A, 

ESR  = ' 1 — 2— : 2 (57) 

A-|  + a2  + A3  + A4 

where  (ESR)1  thru  (ESR)4  are  the  effective  stress  ratios  of  the  four 
triangles.  A-|  thru  A^  are  the  respective  planform  areas  of  the  triangles. 

Now  the  margin  of  safety  MS  is  computed  as  before  by  Eq.  47. 

SHEAR  PANEL 

As  the  name  indicates,  the  shear  panel  is  devised  for  the  purpose  of 
representing  shear  transmitting  elements.  For  example  in  wing  structures 
the  top  and  bottom  skins  can  be  represented  by  membrane  (triangle  and 
quadrilateral)  elements.  If  the  same  elements  are  used  for  spars  and  ribs, 
the  resulting  finite  element  model  grossly  overestimates  the  stiffness 
of  the  structure.  What  this  means  is  that  the  displacements  obtained  by 
this  model  will  be  much  smaller,  or  if  this  model  is  used  for  dynamic 
analysis,  the  frequencies  of  the  structure  will  be  much  higher  and  cannot 
be  matched  with  the  results  obtained  from  ground  vibration  tests.  This 
behavior  is  due  to  the  assumption  of  constant  strain  (stress)  in  the 
membrane  element  formulations.  Most  web  elements  in  box  or  I-beams 
carry  primarily  shear  and  some  normal  stresses.  In  other  words  their 
deformation  is  primarily  due  to  shear  and  not  due  to  normal  stresses. 

The  normal  stresses  in  webs  usually  have  steep  stress  gradients,  and  the 
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assumption  of  constant  stress  (or  strain)  is  not  justified.  To  offset 
this  difficulty,  and  yet  preserve  the  simplicity  of  the  constant  strain 
elements,  a shear  panel  was  formulated  (Reference  8)  with  the  assumption 
that  it  carries  only  shear  stresses.  The  bars  and  other  membrane  elements 
that  surround  the  shear  panel  are  supposed  to  carry  the  normal  stresses. 
Such  a situation  does  not  actually  exist  in  reality,  and  thus  the  shear 
panel  is  an  empirical  element.  However,  the  models  built  on  such  an 
assumption  appear  to  produce  satisfactory  results. 

Until  recently  it  was  a common  practice  in  aircraft  companies  to  model 
wings,  fuselages,  and  empennage  structures  simply  by  bars  and  shear  panels 
to  obtain  primary  load  path  information.  In  such  idealizations  it  was 
a common  practice  to  assign  a third  of  the  cross-sectional  area  as  spar 
and  rib  caps  and  the  remainder  for  the  shear  panels.  It  should  be  pointed 
out  that  every  shear  panel  must  be  surrounded  on  all  four  sides  by  normal 
stress  carrying  elements  such  as  bars  or  membrane  or  bending  elements. 

If  the  natural  model  does  not  contain  such  an  element  on  any  side  of  the 
shear  panel,  a nominal  (or  fictitious)  bar  (post)  must  be  provided. 
Otherwise  the  model  will  have  a singularity. 

The  shear  panel  in  "ANALYZE"  is  constructed  out  of  four  triangles 
with  the  fictitious  node  inside  as  in  the  membrane  quadrilateral  discussed 
earlier.  However,  the  stiffness  matrices  of  the  component  triangles  are 
determined  by  considering  only  the  shear  strain  energy  (Equation  39). 


of  the  element  can  be  sensitive  to  the  reference  axis.  For  rectangular 
elements  the  shear  strain  energy  would  be  the  same  regardless  of  which 
side  is  selected  for  the  reference  axis.  However,  for  quadrilaterals 
the  stiffness  matrix  does  depend  on  the  reference  axis.  The  errors 
produced  by  such  departures  are  usually  not  significant,  but  it  is 
worthwhile  to  make  note  of  the  assumptions  involved.  The  ANALYZE 
program  has  a provision  for  specifying  any  one  of  the  four  sides  of  the 

\ 

quadrilateral  as  the  reference  axis. 

As  in  the  quadrilateral  element  the  shear  stresses  in  all  four 
triangles  are  determined  separately  but  with  respect  to  the  same  reference 
axis.  Of  course,  the  normal  stresses  in  the  shear  panels  have  no 
meaning.  The  margin  of  safety  is  determined  by  a weighted  average  of  the 
effective  stress  ratios(ESR)  as  in  the  quadrilateral.  The  strain  energy 
is  determined  by  considering  only  the  shear  stress  and  strain. 
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4.  ORGANIZATION  OF  THE  PROGRAM 


The  material  presented  in  this  section  is  intended  either  to  help 
introduce  changes  into  the  program  or  to  expand  its  scope  for  the  specific 
needs  of  a researcher  as  the  authors  have  done  in  the  past  ten  years. 

The  steps  outlined  at  the  end  of  Section  2 are  summarized  in  the  flow- 
chart in  Figure  5.  There  are  a total  of  16  boxes  in  the  flow-chart. 

Each  of  these  boxes  generally  involves  one  or  more  subroutines.  The 
subroutines  that  belong  to  each  of  these  boxes  are  identified  first, 
then  the  function  of  each  subroutine  will  be  discussed  in  the  next 
section  with  the  help  of  the  equations  given  in  Sections  2 and  3. 

Box  1 - Input 

Input  in  the  present  version  of  the  "ANALYZE"  program  is  not  in 
subroutine  form.  However,  the  input  statements  are  all  at  the  beginning 
of  the  program,  and  thus  they  can  be  grouped  into  a single  subroutine. 
Alternatively,  one  can  generate  an  input  routine  of  his  own  with 
provisions  like  one  card  per  each  element  and  a card  for  each  node  etc. 

For  example,  it  is  relatively  easy  to  write  a subroutine  with  NASTRAN  type 
input.  The  description  of  the  various  arrays  (See  input  instructions) 
and  their  dimension  requirements  given  in  Appendix  A can  be  quite  helpful 
in  writing  such  an  input  routine. 

Box  2 - Map  Stiffness  Matrix 

This  step  involves  a single  subroutine  called  "POP".  The  purpose  of 
this  routine  is  simply  to  estimate  the  storage  requirements  of  the  stiffness 
matrix  and  to  map  its  profile.  The  stiffness  matrix  is  stored  in  a single 
array  called  SK.  The  elements  of  the  matrix  are  stored  columnwise 
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Figure  5.  Flow  Chart  for  Program  "ANALYZE 


starting  from  the  first  non-zero  element  in  the  column  to  the  diagonal 
element.  Since  the  matrix  is  symmetric,  only  the  upper  triangle  is 
stored. 

Box  3 - Element  Stiffness 

There  are  four  elements  in  the  program.  All  of  them  require  the 
subroutine  "COORD".  In  addition  all  the  plate  elements  require  the 
routine  "ELSTIC".  The  remaining  subroutines  are  listed  separately  for 
each  element. 

i.  Bar  (Rod)  Element: 

The  bar  element  is  shown  in  Figure  2a  with  the  local  coordinate 
system  and  degrees  of  freedom.  This  element  requires  the  subroutine 
"ELSTIF"  which  generates  the  bar  stiffness  matrix  in  the  local  coordinate 
system  and  also  transforms  it  to  the  global  coordinate  system. 

ii.  Triangular  Membrane  Element: 

The  element  and  its  local  coordinate  system  are  shown  in  Figure  2b. 
The  subroutine  "PLSTIF"  is  the  only  other  routine  required  by  this  element 
It  generates  the  stiffness  matrix  of  the  triangle  in  the  local  coordinate 
system. 

iii.  Quadrilateral  Membrane  Element  and  Shear  Panel 

The  elements  and  their  local  coordinate  system  are  shown  in  Figure  2c 
The  subroutines  "QDRLTL" , "PLSTIF",  "SUM",  "CONDNS",  "CHANGE"  and  "CRAMER" 
are  the  additional  routines  required  by  these  elements.  Together  these 
subroutines  generate  the  stiffness  matrix  of  either  the  quadrilateral 
membrane  or  shear  panel.  The  routine  "QDRLTL"  calls  "PLSTIF",  "SUM" 
and  "CONDNS".  The  routine  "PLSTIF"  calls  "CRAMER".  Similarly  "CONDNS" 
calls  "CHANGE". 
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Box  4 - Transform 


This  step  involves  a single  subroutine  called  "TRNSFM".  It 
transforms  the  stiffness  matrices  of  the  triangles,  quadrilaterals,  and 
shear  panels  from  the  local  to  the  global  coordinate  system. 

Box  5 - Assemble 

"ASEMBL"  is  the  only  subroutine  used  in  this  step.  Its  purpose  is 
to  add  the  element  stiffness  matrices  to  the  total  stiffness  matrix  of 
the  structure.  The  steps  3 thru  5 form  a lo^o  in  which  all  the  element 
stiffness  matrices  are  computed  and  assembled  into  the  total  stiffness 
matrix. 

Box  6 - Boundaries 

The  routine  called  "B0UND2"  eliminates  the  rows  and  columns  of  the 
stiffness  matrix  corresponding  to  the  support  degrees  of  freedom  of  the 
structure.  In  addition  it  also  condenses  the  stiffness  matrix. 

Box  7 - Reduce  Force 

This  step  involves  a routine  called  "REDUCE".  It  eliminates  the  rows 
of  the  force  matrix  corresponding  to  the  support  degrees  of  freedom. 

Box  8 - Solution  of  the  Force  Deflection  Equations 

The  routine  "GAUSS"  solves  the  load  deflection  equations  by  Gaussian 
elimination.  A large  percentage  of  the  analysis  time  (80  to  90%)  is 
spent  in  this  routine,  and  its  efficiency  is  extremely  important  in  reducing 
the  costs  of  the  analysis.  At  the  end  of  this  step  the  displacements  of 
the  structure  are  available  in  condensed  form  (excluding  boundary  degrees 
of  freedom)  in  the  global  coordinate  system. 
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Box  9 - Total  Energy 


The  total  energy  of  the  structure  is  computed  by 
W = \ R*  r 

The  strain  energy  of  the  structure  (U)  is  computed  by  adding  the  strain 
energies  of  the  elements  in  step  14  (Box  14).  A comparison  of  W and  U 
provides  an  equilibrium  check. 

Boxes  10  and  11  - Restore  Forces  and  Displacements 

These  two  boxes  use  the  same  routine  called  "RESTOR".  The  purpose  of 
this  routine  is  to  restore  the  force  and  deflection  matrices  to  their 
original  dimension  to  include  the  boundary  degrees  of  freedom.  Its 
purpose  is  essentially  opposite  to  that  of  the  routine  "REDUCE"  in  Box  7. 

Box  12  - Element  Displacements 

The  routine  "COORD"  and  "ELFORC"  facilitate  extraction  and  transformation 
of  the  element  displacements  from  the  global  to  the  local  coordinate  system. 
Box  13  - Element  Forces 

This  step  is  not  in  all  versions  of  "ANALYZE".  Element  forces  are  not 
necessary  to  compute  stresses.  However,  this  step  can  be  restored  if  the 
element  shear  flows  and  other  force  information  are  necessary. 

Box  14  - Element  Stresses 

The  details  of  this  step  depend  on  the  type  of  element, 
i.  Bar  (Rod)  Element: 

The  stress  in  this  element  is  computed  in  the  program  itself.  No 
additional  routines  are  involved.  At  the  same  time  the  element  strain 
energy  is  also  computed. 
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ii.  Triangular  Membrane  Element: 

The  subroutines  "STRESS"  and  "CRAMER"  are  involved  in  this  step.  The 
routine  "STRESS"  calls  "CRAMER".  The  purpose  of  this  routine  is  to  calculate 
stresses  in  the  triangular  element.  In  addition  this  routine  calculates 
strain  energy  and  the  effective  stress  in  the  element  (See  Equations 
44  and  45). 

iii.  Quadrilateral  Membrane  and  Shear  Panel 

This  step  involves  routines  "ELSTIC",  "QDRLTL" , "PLSTIF",  "SUM", 
"CONDNS",  "CRAMER",  "QLSTRS"  and  "STRESS".  It  should  be  noted  that  the 
routine  "QDRLTL"  calls  "PLSTIF",  "SUM  and  "CONDNS".  "PLSTIF"  in  turn 
calls  "CRAMER". 

Box  15  - Stress  Output 

The  instructions  for  the  output  of  the  table  of  stresses  are  in  the 
main  program.  No  subroutine  is  used  for  the  output  itself.  The  steps 
12  thru  15  form  a loop  in  which  the  stress  information  for  all  the  elements 
is  computed  and  printed  in  a table.  This  is  one  of  the  two  main  tables 
of  output  of  this  program.  Explanation  of  this  table  is  given  in  the 
section  on  output  (Section  7). 

Box  16  - Displacement  Output 

This  step  involves  a single  subroutine  called  "PRNTDR".  This  routine 
prints  out  the  second  important  table  of  output  which  contains  information 
about  the  nodes.  This  information  includes  the  coordinates  of  the  nodes, 
applied  forces  and  the  calculated  displacements  for  each  node.  The 
detailed  explanation  of  this  table  is  given  in  the  section  on  output 
(Section  7). 


In  addition  to  the  above  16  steps  there  are  instructions  for  weight 
computations  and  other  details,  and  their  purpose  can  be  identified  from 
the  program.  There  are  very  few  comment  cards  in  the  main  body  of  the 
program  and  this  omission  is  by  design  in  order  to  avoid  continuous 
updating.  The  user  can  incorporate  his  own  cogent  cards  with  the 
help  of  the  explanation  given  in  this  section. 
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5.  DESCRIPTION  OF  THE  SUBROUTINES 


"ANALYZE"  consists  of  the  main  program  and  21  Subroutines,  The  main 
program  has  260  cards.  The  length  of  the  Subroutines  varies  from  15  to 
62  cards.  The  total  length  of  the  program  is  under  1000  cards.  A list 

of  the  Subroutines,  the  number  of  Cards  in  each  Subroutine  and  other 

details  are  given  in  Table  1.  The  flow  chart,  Fig.  5,  and  the  explanation 
in  the  previous  section  give  details  of  the  main  program.  The  description 
of  the  Subroutines  is  given  in  the  remainder  of  this  section. 

Subroutine  "POP" 

The  purpose  of  Subroutine  "POP"  is  to  estimate  the  storage  requirements 

of  the  stiffness  matrix  before  actually  determining  it.  This  information 

can  be  generated  from  the  element  connections  with  the  nodes.  For  example, 
if  an  element  connects  4 nodes,  and  if  each  node  has  3 degrees  of  freedom 
in  the  global  coordinate  system,  then  the  stiffness  matrix  of  the  element 
would  be  of  dimension  12  x 12.  This  matrix  can  be  partitioned  four  ways, 
in  both  row  and  column  directions  as  shown  in  Fig.  6.  The  location  of 
these  sixteen  submatrices  in  the  total  stiffness  matrix  can  be  determined 
by  the  address  of  the  nodes  to  which  the  element  is  connected.  If  the 
element  is  connected  to  the  nodes  MA,  MB,  MC,and  MD,  then  the  addresses  of 
the  element  submatrices  in  the  total  stiffness  matrix  are  shown  in  Fig.  6. 

If  all  the  elements  are  connected  to  all  the  nodes,  then  the  stiffness 
matrix  of  the  structure  will  be  fully  populated.  The  non-zero  elements 
in  the  matrix  are  considered  as  population.  Since  most  of  the  elements 
connect  only  a few  nodes,  the  stiffness  matrices  are  usually  sparsely 
populated.  Determining  the  profile  of  the  stiffness  matrix  population 
is  the  essential  function  of  the  routine  "POP". 
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Total  Stiffness  Matrix 
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The  distribution  of  the  nonzero  elements  is  dependent  upon  the 

way  the  nodes  of  the  finite  element  model  are  numbered.  Because  of  the 

symmetry  of  the  stiffness  matrix,  only  the  lower  or  upper  triangular 

matrix  is  considered.  For  the  purpose  of  this  discussion  definitions  of 

the  following  terms  are  in  order.  The  gross  population  (Pgross)  of  the 

stiffness  matrix  is  defined  as  the  total  number  of  elements  in  the 

upper  triangle  of  the  matrix.  The  net  population  (Pnet)  is  the  total 

number  of  non-zero  elements  in  the  upper  triangle.  Zeros  resulting  from 

transformations  are  not  excluded  from  the  net  population.  The  apparent 

population  (P  is  the  actual  number  of  elements  considered  as 

apparent 

nonzeros  by  a given  solution  scheme.  From  these  definitions 


net 

For  a given  structure  P 


apparent 


gross 


„„  and  P . are  invariant  and  are  given  by 
gross  net  3 J 


(59) 


gross 


N (N  + 1 ) 
2 


(60) 


and 


m « 

pnet  = n (n  + 1 ) (number  of  nodes)  + E n Lk.j  (k^.-l )] 

2 l—l  M 


(61) 


where  N is  the  total  number  of  degrees  of  freedom  of  the  structure,  n is 
the  number  of  degrees  of  freedom  of  each  node  (all  the  nodes  are  assumed 
to  have  the  same  number  of  degrees  of  freedom;  when  this  is  not  true  the 
necessary  modification  is  simple),  k.  is  the  number  of  nodes  to  which  the 
ith  element  is  connected,  and  m is  the  number  of  elements  in  the  structure. 
The  quantity  NR  is  given  by 


NR  = Z (b.  - 1) 
1=1  1 


(62) 
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INTERSECTING  PLATES 


FIGURE  7:  DISTRIBUTION  OF  NONZERO  ELEMENTS  IN  THE  STIFFNESS  MATRIX 


Subroutine  "ELSTIC 


This  routine  generates  the  3 x 3 elastic  constants  matrix  for  a 
given  material  (see  Eq.  3). 

Subroutine  “COORD" 

This  routine  establishes  the  local  coordinate  system  for  all  the 
elements  and  also  determines  the  nodal  coordinates  in  the  local  system. 
It  generates  the  direction  cosine  matrix  which  will  be  used  to  transform 
the  element  stiffness  matrices  to  the  global  coordinate  system  (see 
Eqs.  13  and  16). 

i.  Bar  Element 

The  local  coordinate  system  of  the  bar  element  is  established  by 
drawing  a line  between  the  two  nodes  MA  and  MB  (see  Fig.  2)  connecting 
the  bar.  The  direction  cosines  are  determined  by 

XComp  = XMA  " XMB 


Y = Y - Y 
Comp  ’MA  ’MB 


ZComp  ~ ZMA  " ZMB 


L = (XZ  + YZ  + ZZ  )^ZZ 
'Comp  Comp  Comp' 


(64) 


(65) 


Y Z 

Comp  _ Comp 

L nl  ” L 


(66) 


where  X^,  Y^  and  Z^  are  the  three  coordinates  of  the  node  MA  in  global 
coordinate  system.  The  direction  cosines  , m^and  n^  become  the  first 
row  of  the  3x3  matrix  A. 
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11.  Triangular  Membrane  Element 

The  local  coordinate  system  of  the  triangular  membrane  element  is 
established  by  assigning  the  local  x-axis  to  the  line  joining 
nodes  MA  and  MB.  The  direction  cosines  of  this  line  are  determined  as 
in  the  case  of  the  bar  element.  The  plane  of  the  plate  is  established  by 
two  unit  vectors  in  the  directions  of  the  lines  joining  nodes  MA-MB  and 

A /\ 

MA-MC.  If  a and  b are  these  two  unit  vectors,  then  the  normal  to  the 


plane  is  obtained  by 


a x b = c 


Since  a and  b are  not  orthogonal  vectors,  c is  not  a unit  vector. 
The  unit  vector  in  this  direction  is  given  by 


A 

The  local  z-axis  is  in  the  direction  of  the  unit  vector  c.  Now 
the  local  y-axis  is  established  by 

c x a = d (69) 

The  direction  cosines  of  x and  y become  the  first  two  rows  of  matrix  A. 

iii.  Quadrilateral  Membrane  and  Shear  Panel 

The  local  coordinate  system  of  the  quadrilateral  membrane  and  the 
shear  panel  are  established  by  a procedure  similar  to  that  of  the  triangle. 
The  plane  of  the  triangle  connecting  the  three  nodes  MA,  MB,  and  MC  becomes 
the  reference  plane.  Any  warping  in  the  quadrilaterals  and  shear  panels  is 
ignored.  If  there  is  too  much  warping  in  the  quadrilaterals,  it  is  better 
to  divide  them  into  two  or  more  triangles  or  reduce  the  mesh  size.  In  the 
case  of  excessively  warped  shear  panels,  the  size  of  the  grid  must  be 


40 


reduced.  "ANALYZE"  does  not  have  a provision  for  determining  the  warp 
and  the  consequent  kick  forces. 

The  node  MA  of  the  element  becomes  the  origin  of  the  element  local 
coordinate  system  and  the  coordinates  of  the  remaining  nodes  are  determined 
by  expressions  similar  to  the  following: 

X3  = ^ XMC  " W£1  + (YMC  " Wml  + ' Wnl 

y3  = ^ XMC  ‘ W£2  + ^YMC  " YMA^m2  + ^ ZMC  ' Wn2 

This  subroutine  also  determines  the  coordinates  of  the  fictitious  node 
needed  to  break  the  quadrilateral  and  shear  panels  into  four  triangles. 

This  interior  node  is  established  by 

X1  + x2  + x3  + x4 

x5  4 

(70) 

= VVV_y4 
y5  4 


where  x-j , x?  xg  and  y-j , y^  yg  are  the  coordinates  of  the  five  nodes 

(including  the  fictitious  interior  node)  of  the  quadrilaterals  and  shear 
panels  in  the  local  coordinate  system. 


Subroutine  "ELSTIF" 

This  subroutine  determines  the  stiffness  matrix  of  the  bar  by  Eq.  22. 
It  also  transforms  the  bar  stiffness  matrix  to  the  global  coordinate  system 
by 

Ki  = af  ^ a5  (71) 
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Subroutines  "PLSTIF"  and  “CRAMER 


The  routine  "PLSTIF"  determines  the  element  stiffness  matrix  of  the 
triangle  in  the  local  coordinate  system.  This  is  also  the  basic  routine 
for  determining  the  stiffness  matrices  of  the  four  triangles  of  the 
quadrilateral  and  the  shear  panel. 


"PLSTIF"  first  calls  the  routine  "CRAMER",  which  determines  the 
inverse  of  the  matrix  X by  Cramer's  rule.  The  matrix  X is  given  by  Eq.  34. 
The  determinant  of  X represents  twice  the  area  of  the  triangle. 

Then  the  "PLSTIF"  subroutine  determines  the  element  stiffness  matrix 
by  Eq.  40.  In  determining  the  matrices  and  e^,  it  takes  advantage 
of  the  fact  that  the  columns  of  Z~^  (see  Eq.  33)  represent  unit  displacement 
modes  (see  explanation  under  Eq.  34). 


In  computing  the  stiffness  matrices  of  the  triangles  of  the  shear 
panels,  "PLSTIF"  considers  only  the  shear  strain  energy.  For  example,  in 
such  a case,  Eq.  40  becomes 


* = i 


(i)  0) 

exyGexy 


(6)  (1) 

exyGexy 


(1)  (2) 

exyGexy“* 


(6)  (2) 
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xy  xy 
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Subroutine  "QDRLTL" 

This  subroutine  simply  manages  the  routines  "PLSTIF",  "SUM",  and 
"CONDNS"  in  computing  the  stiffness  matrix  of  the  quadrilateral  membrane 
and  shear  panel.  This  routine  also  makes  provision  for  assigning  different 
sides  as  reference  axis  for  the  shear  panels. 


Subroutine  "SUM" 

This  subroutine  adds  the  four  triangle  stiffness  matrices  computed 
by  "PLSTIF"  to  produce  a 10  x 10  stiffness  matrix  (including  two  degrees 
of  freedom  for  the  interior  node)  for  the  quadrilateral  or  shear  panel. 

Subroutine  "CONDNS" 

This  routine  condenses  the  10  x 10  quadrilateral  or  shear  panel 
stiffness  matrix  to  an  8 x 8 matrix.  The  condensation  is  done  by  using 
Eq.  56. 

Subroutine  "CHANGE" 

This  routine  interchanges  the  rows  and  columns  of  the  quadrilateral 
(or  shear  panel)  stiffness  matrix  so  that  the  element  degrees  of  freedom 
are  in  ascending  order  before  addition  to  the  structure  stiffness  matrix. 
This  step  is  necessary  because  the  routine  "ASEMBL"  assumes  that  the 
element  degrees  of  freedom  are  In  ascending  order. 

Subroutine  "TRNSFM" 

This  routine  transforms  the  plate  element  stiffness  matrices  from  the 
local  to  the  global  coordinate  system  by  (see  Eq.  16) 

Ki  = aj  ^ ?1  (73) 

where  is  the  transformed  element  stiffness  matrix  of  the  ith  element. 

Subroutine  "ASEMBL" 

This  routine  adds  the  element  stiffness  matrices  to  the  total  stiff- 
ness matrix. 
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For  an  explanation  of  the  rules  of  this  addition  see  the  description  of 
subroutine  "POP".  It  should  be  noted  that  only  the  upper  half  of  the 
stiffness  matrix  is  stored.  This  storage  is  columnwise  starting  with  the 
first  non-zero  element  above  the  diagonal. 


Subroutine  "PRINTK" 

The  purpose  of  this  routine  is  to  print  the  stiffness  matrix  (if 
desired)  rowwise  starting  with  the  first  non-zero  element  and  proceeding 
to  the  diagonal. 

Subroutine  "B0UND2" 

This  routine  eliminates  the  rows  and  columns  corresponding  to  the 
constrained  degrees  of  freedom  and  condenses  the  stiffness  matrix. 

Subroutine  "REDUCE" 

This  routine  eliminates  the  rows  of  the  applied  force  matrix 
corresponding  to  the  constrained  degrees  of  freedom.  It  is  assumed  that 
i ch  column  of  the  force  matrix  represents  an  independent  load  condition. 

Suoroutine  "GAUSS" 

"GAUSS"  solves  the  load  deflection  equations  (Eq.  17)  by  Gaussian 
elimination.  The  first  step  of  the  solution  is  the  decomposition  of  the 
stiffness  matrix  by  Eq.  18.  The  next  two  steps  represent  forward  and  back 
substitution  using  Eqs.  19  and  20  respectively.  For  the  solution  of 
additional  load  vectors  only  the  steps  FBS  have  to  be  repeated.  If  "GAUSS" 
is  entered  with  any  value  other  than  0 for  the  parameter  NDCOMP,  only  the 
last  two  steps  will  be  executed.  The  matrices  L and  D are  stored  in  place 
of  the  original  stiffness  matrix. 
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Subroutine  "RESTOR 


This  routine  restores  the  displacement  or  force  matrix  to  full  size 
by  assigning  zero  values  to  boundary  degrees  of  freedom. 

Subroutine  "ELFORC" 

This  routine  extracts  the  element  displacements  from  the  global 
coordinate  system  and  transforms  them  to  the  local  coordinate  system  by 
Eq.  13. 

Subroutine  "STRESS" 

The  purpose  of  the  "STRESS"  routine  is  to  compute  strains  and  stresses 
in  the  triangular  element.  It  first  calls  the  routine  "CRAMER"  which 
computes  X ^ (Eq.  34)  by  Cramer's  rule.  The  strains  in  the  element  are  then 
calculated  by  Eqs.  30  and  35  thru  37.  The  stresses  in  the  element  are 
computed  by  Eq.  2.  Also  it  computes  the  strain  energy  and  the  effective 
stress  in  the  element  by  Eqs.  1 and  45  respectively. 

Subroutine  "QLSTRS" 

This  routine  prepares  the  data  for  computing  stresses  in  the  four 
triangles  of  the  quadrilateral  or  shear  panel  elements.  First  it  determines 
the  interior  node  displacements  from  the  corner  node  displacements  using 
Eq.  54.  Then  it  calls  subroutine  "STRESS"  to  compute  the  stresses  in  the 
four  triangles.  It  adds  the  strain  energy  of  the  four  triangles  to  obtain 
the  total  strain  energy.  It  identifies  the  triangle  with  the  largest 
effective  stress  and  normalizes  the  effective  stress  of  the  three  remaining 
triangles  with  respect  to  this  largest  value. 

Subroutine  "PRNTDR" 

This  subroutine  prints  out  the  table  of  node  information.  This  includes 
the  node  number,  its  coordinates,  applied  forces,  and  the  displacements. 
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NAME 

NUMBER  OF  CARDS 

CALLED  FROM 

ANALYZE 

315 

Main  Program 

POP 

62 

ANALYZE 

ELSTIC 

15 

ANALYZE 

COORD 

44 

ANALYZE 

ELSTIF 

21 

ANALYZE 

PLSTIF 

46 

ANALYZE,  QDRLTL 

CRAMER 

19 

PLSTIF,  STRESS 

QDRLTL 

32 

ANALYZE 

SUM 

23 

QDRLTL,  QLSTRS 

CONDNS 

36 

QDRLTL,  QLSTRS 

CHANGE 

25 

CONDNS 

TRNSFM 

36 

ANALYZE 

ASEMBL 

41 

ANALYZE 

PRINTK 

15 

ANALYZE 

BOUND2 

' 35 

ANALYZE 

REDUCE 

18 

ANALYZE 

GAUSS 

57 

ANALYZE 

RESTOR 

28 

ANALYZE 

ELFORC 

22 

ANALYZE 

STRESS 

33 

ANALYZE,  QLSTRS 

QLSTRS 

65 

ANALYZE 

PRNTDR 

39 

ANALYZE 

TOTAL 

1027 

Table  1:  Program  Description 
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6.  INPUT  INSTRUCTIONS 


Input  for  the  programs  is  divided  into  a number  of  card  sets. 
Each  card  set  will  consist  of  one  or  more  cards.  Only  three 
Formats  are  used  for  input.  An  integer  Format  (1415),  a floating 
point  Format  (6F10.0)  and  a mixed  Format  3(F10. 0,215).  The  first 
four  card  sets  will  each  have  one  card  regardless  of  the  size  of  the 
problem.  The  number  of  cards  required  for  the  remaining  card  sets 
depends  on  the  problem  size.  The  first  card  set  indicates  the 
number  of  problems  (structures)  to  be  analyzed.  If  this  number  is 
more  than  one,  the  program  assumes  that  the  remaining  card  sets  will 
be  supplied  for  each  problem  one  after  the  other.  The  next  card  set 
is  for  the  title  of  the  problem.  Card  set  three  defines  the  basic 
parameters  like  the  number  of  elements,  nodes  etc.  And  set  4 defines 
the  properties  of  a reference  material.  This  material  can  be  any  one 
of  the  materials  used.  The  remaining  card  sets  define  material 
properties  (5  and  6),  type  of  elements  (7),  element  connections 
(8,  9,  10,  11),  sizes  o^  the  elements  (12),  element-material 
identification  (13),  node  coordinates  (14),  boundaries  (15)  and 
loading  information  (16  and  17). 
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INPUT  INSTRUCTION  DETAILS 


CARD  SET 
(FORMAT) 

1 

(1415) 

2 

(8A10) 

3 

(1415) 


4 

(6F10.5) 


PARAMETER 

DESCRIPTION 

NSTR 

Number  of  data  sets 

TITLE 

An  alphanumeric  description  of  the 
problem  to  be  solved. 

MEMBS 

Number  of  elements 

JOINTS 

Number  of  nodes 

NBNDRY 

Number  of  restrained  degrees 
of  freedom 

LOADS 

MM 

Number  of  loading  conditions 

m\=2  Two  dimensional  problem 
nL=3  Three  dimensional  problem 

NR 

Variable  used  only  for  calculating 

the  net  population  of  the  stiffness 
matrix.  It  has  no  other  role  in 
the  program.  See  Section  5. 


INCHES 

TNrHF<.f=l  Coordinate  data  is  in  inches 
LH/1  Coordinate  data  is  in  feet 

KIPS 

r=l  Applied  forces  are  in  kips 

Kil/i  Applied  forces  are  in  pounds 

NMAT 

Number  of  materials 

MSSTRS 

r=0  Margin  of  safety  calculated  from 
M__TR_  default  allowable  stresses. 

Margin  of  safety  calculated  from 
input  allowable  stresses. 

EEE 

YOUNG'S  modulus/106  of  one  of  the 
elements  in  psi. 

PMU 

POISSON'S  ratio  of  one  of  the 
elements. 

RHO 

Density  of  one  of  the  elements  in 
lbs/in3. 
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CARD  SET 
(FORMAT) 

PARAMETER 

DESCRIPTION 

IF  MSSTRS 

= 0,  skip  CARD  SET  5. 

5 

(6F10.5) 

ALSTRS(I) 

1=1 , . . . ,3*NMAT 

3 

Allowable  st.resses/10  in  tension, 
compression  and  shear  for  the  1™ 
material . 

IF  NMAT  f 
IF  NMAT  = 

1 , CARD  SET  4 parameters 
1,  skip  CARD  SET  6. 

can  be  for  any  of  the  materials. 

6 

(6F10.5) 

YOUNGM(I) 

YOUNG'S  modulus/106  for  the  Ith 
material  in  psi. 

P0 I SON ( I ) 

POISSON'S  ratio  for  the  It,n  material. 

RH01(I) 

1=1,..., NMAT 

Density  for  the  Ith  material  in 
lbs/in3. 

7 

(1415) 

NNODES ( I ) , 1=1 , . . . ,MEMBS  Element  Type 

r=2  BAR 

NNODE(l)  TRIANGLE 

nnuutu;  >=4  QUADRILATERAL  MEMBRANE 

-5  SHEAR  PANEL 

8 

(1415) 

MA ( I ) , 1=1 ,. . . ,MEMBS 

First  node  number  of  each  element. 

9 

0415) 

MB  (.1 ) ,1=1 ....  ,MEMB5 

Second  node  number  of  each  element. 

10 

(1415) 

MC( I ) , 1=1 .... ,MEMBS 

Third  node  number  of  each  element. 

11 

(1415) 

MD( I ) , 1=1 , . . . ,MEMBS 

Fourth  node  number  of  each  element. 

NOTE:  For  bars  leave  MC(I)  and  MD(I)  blank.  For  triangles  leave  MD(I)  blank. 

For  each  element  let  MA(I)  be  the  lowest  node  number  and  MB(I)  be  the 
next  lowest.  For  Quadrilaterals  and  Shear  Panels,  MC(I)  and  MD(I)  are 
determined  by  continuing  in  the  direction  defined  by  MA(I)  and  MB ( I ) . 

12 

(6F10.5) 

TH ( I ) ,1=1 .... ,MEMBS 

Thickness  of  each  element. 

For  a bar  thickness  is  cross-sectional 
area. 

IF  NMAT  =1,  skip  CARD  SET  13. 

13 

(1415) 

MYOUNG(I),  1=1 MEMBS 

Material  number  of  each  element. 
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7.  OUTPUT  DESCRIPTION 

The  primary  output  of  the  program  ANALYZE  consists  of  two 
tables  (items  6 and  8 of  the  output  description  details).  The 
first  table  gives  element  information  and  the  second  table  gives 
information  about  the  nodes.  The  element  information  includes 
member  number,  thickness  (cross-sectional  area  of  the  bars), 
planform  area  (length  of  a bar),  element  type,  stress  information, 
strain  energy,  and  margin  of  safety.  The  information  about  the 
nodes  includes  node  (joint)  number,  node  coordinates,  applied  forces, 
and  the  resulting  displacements.  In  addition  to  these  two  tables 
output  3a  (coming  from  subroutine  POP)  gives  important  information 
about  the  population  distribution  of  the  stiffness  matrix.  The 
value  of  the  apparent  population  is  crucial  in  determining  the 
dimension  of  the  stiffness  matrix  (SK).  This  dimension  must  be  at 
least  as  big  as  or  bigger  than  this  value. 

Item  7 gives  information  about  the  total  strain  energy  (U) 
and  the  work  of  the  external  forces  (W)  for  the  structure.  This 
information  can  be  very  useful  for  an  equilibrium  check. 

Item  4 gives  the  weight  of  the  structure.  The  remaining 
information  is  not  really  very  Important  to  the  user. 


OUTPUT  DESCRIPTION  DETAILS 


Output  for  Program  ANALYZE  consists  of  the  following: 

1)  Untitled  echo  of  all  input  data  except  boundaries  and  applied 
loads. 

2)  Boundary  data,  i.e.  contents  of  array  IBND  (CARD  SET  15) 

3)  Output  from  Subroutine  POP  concerning  the  distribution  of  elements 
in  the  stiffness  matrix.  This  information  is  generated  before  the 
stiffness  matrix  of  the  structure  is  assembled. 

(a)  Gross  Population  - total  number  of  elements  in  the  upper 
triangle  of  the  matrix. 

Net  Population  = actual  population  of  possible  non-zero  elements 
in  the  upper  triangle  of  the  stiffness  matrix.  This  number  would 
be  correct  only  if  NR  is  correct  in  CARD  SET  3. 

Apparent  Population  = actual  number  of  elements  considered  as 
non-zero  by  a given  solution  scheme.  Thus  the  apparent 
population  represents  the  number  of  storage  locations  required 
for  the  stiffness  matrix. 

(b)  Starting  Row  Numbers  for  each  column  - the  number  of  the  row 
where  the  first  non-zero  element  occurs  in  each  column. 

(c)  Number  of  Diagonal  Elements  in  Single  Array  Stiffness  Matrix. 

For  each  Column  I the  actual  number  of  elements,  ID( I ) , in  the 
upper  triangular  matrix  up  to  and  including  that  column,  i.e. 

nm-JUta  - l b, 

* j=l  J 

where  bj  is  the  row  number  given  for  Column  I in  (b).  Thus 
for  the  last  column,  ILAST, 

ID(ILAST)  = Apparent  Population 

4)  Weight  of  the  structure 

5)  Boundary  conditions  applied  to  the  stiffness  matrix  alter  the  arrays 
defined  by  (b)  and  (c)  above,  and  thus  they  are  reprinted. 

6)  Output  for  each  element  after  analysis. 

(a)  MEMB  - Element  Number 

(b)  THICK  - Thickness  of  the  element.  For  a bar  thickness  is 

cross-sectional  area. 
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(c)  AREA  - Area  of  the  element.  For  a bar  area  is  length. 

(d)  TYPE  - Type  is  a composite  number  which  describes  the  element 
type  and  material  number.  Type  is  defined  as 

TYPE  = NNODES ( I )xlO  + MYOUNG(I). 

See  CARD  SETS  7 and  13. 

Note;  If  the  number  of  materials  is  greater  than  10,  TYPE  is 
meaningless.  If  the  number  of  materials  is  1,  MY0UNG(I)=1  for 
all  I. 

(e)  MA,  MB,  MC,  MD  - defined  in  CARD  SETS  8,9,10,and  11. 

(f)  SIGMA-X  (a  ) , SIGMA-Y  (a),  SIGMA-XY  (a  ). 

x V xy 

Stresses  in  the  x-y  local  coordinates  of  the  element. 

EFSTR-1 , EFSTR-2,  EFSTR-3,  EFSTR-4  - Effective  stresses 
in  the  element  determined  by  the  Von  Mises  Criterion. 

The  stress  output  varies  per  element  type. 

(i)  BAR  SIGMA-X  only 

(ii)  TRIANGLE  SIGMA-X,  SIGMA-Y,  SIGMA-XY,  EFSTR-1 

(iii ) QUADRILATERAL  MEMBRANE 

The  Quadrilateral  membrane  element  is  divided  into 
4 triangles  for  analysis.  SIGMA-X,  SIGMA-Y,  SIGMA-XY  are 
for  that  triangle  with  the  maximum  effective  stress.  This 
maximum  effective  stress  is  given  as  EFSTR-i  for  some  i, 
i=l,...,4.  Then  EFSTR-j,  j^i,  are  defined  as  the  ratio  of 
the  effective  stress  for  triangle  j to  the  maximum  effective 
stress. 

(iv)  SHEAR  PANEL 

The  Shear  Panel  is  also  divided  into  4 triangles  for  analysis. 

SIGMA-XY  ( t ) is  for  that  triangle  with  the  maximum  effective 
xy 

stress.  Then  EFSTR-i,  i=l  »• . ..,4  are  as  defined  in  (iii). 

(g)  ENERGY  - Total  strain  energy  in  the  element. 

(h)  MS  - Margin  of  Safety  for  the  element. 

NOTE;  If  the  number  of  loading  conditions  is  greater  than  1,  output  (f) 
and  (h)  are  given  continuously  for  each  load  case. 

7)  The  total  strain  energy  (U)  of  the  structure  and  the  work  (W)  of  the 
external  forces  for  each  loading  condition. 
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8)  Output  for  each  node  after  analysis. 


(a)  JOINT  - Node  Number 

(b)  X,  Y,  Z - x,  y,  and  z coordinate  of  the  node 

(c)  FORCE-X,  FORCE-Y,  FORCE-Z  - applied  forces  in  the  x,  y,  and 
z direction. 

(d)  DISPL-X,  DISPL-Y,  DISPL-Z  - Displacements  in  the  x,  y, and  z 
direction. 

NOTE:  If  the  number  of  loading  conditions  is  greater  than  1,  output 
(c)  and  (d)  are  given  continuously  for  each  load  case. 


8.  SAMPLE  PROBLEM 


A three  spar  wing  shown  in  Figure  8 is  idealized  by  membrane 
quadrilaterals,  shear  panels, and  bars  (axial  force  members).  The  top 
and  bottom  skins  are  represented  by  membrane  quadrilaterals  and  the 
spars  and  ribs  by  shear  panels.  In  addition,  the  top  and  bottom  nodes 
are  connected  by  bar  elements  or  posts.  The  root  section  of  the  wing 
is  assumed  to  be  fixed.  The  finite  element  model  of  the  wing  box  is 
shown  in  Figure  8. 

The  wing  is  analyzed  for  two  independent  loads.  These  loading 
conditions  are  generated  by  simplified  pressure  distributions  representative 
of  a subsonic,  forward-center-of-pressure  loading, and  a supersonic  near- 
uniform-pressure loading.  The  material  properties  (aluminum)  are  given 
in  Figure  9. 

The  results  of  the  analysis  are  given  in  Appendix  D.  They  are  given 
in  two  tables.  The  first  table  gives  element  information.  This  includes 
sizes,  connections,  stresses,  energies,  and  margins  of  safety.  The  second 
table  gives  information  about  nodes.  This  includes  coordinates,  applied 
loads, and  displacements.  For  a detailed  description  of  the  complete  output, 
see  Section  7. 

The  wing  was  also  analyzed  by  the  NASTRAN  program.  Table  2 compares 
ANALYZE  and  NASTRAN  z-displacements. 
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Figure  8.  Aerodynamic  Planform  and  Primary  Structural  Arrangement  of  Wing 
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APPENDIX  A:  ESTIMATION  OF  CORE  REQUIREMENTS 


The  purpose  of  this  appendix  is  to  aid  in  the  approximate 
estimation  of  core  requirements  for  the  program.  These  change  with 
the  problem  size.  For  example  with  100kg  (120K  bytes  or  30K  decimal) 
it  is  possible  to  solve  problems  of  the  size  250  to  300  degrees  of 
freedom  assuming  that  the  nodes  are  numbered  with  reasonable  care  for 
an  optimum  stiffness  profile  (See  the  discussion  under  subroutine  "POP" 
in  Section  5).  The  dimensional  requirements  of  various  arrays  are  explained 
by  comment  cards  at  the  beginning  of  the  program.  However,  this  section 
reiterates  the  importance  of  adjusting  the  dimensions  of  some  important 
arrays. 

The  arrays  can  be  grouped  into  nine  types.  The  number  of  elements, 
degrees  of  freedom,  loading  conditions,  the  number  of  boundaries  and 
the  number  of  materials  are  some  of  the  variables  that  affect  the  size 
of  the  arrays.  The  arrays  must  be  dimensioned  at  least  as  big  or 
bigger  than  the  number  of  these  variables  in  the  problem.  The  arrays 
with  fixed  sizes  (not  affected  by  problem  size)  are  dimensioned  first. 

The  total  core  requirement  of  these  arrays  is  relatively  small.  Next, 
the  arrays  that  depend  on  the  number  of  materials  are  dimensioned. 

The  third  group  consists  of  a single  array  IBND  which  is  dimensioned 
according  to  the  number  of  boundary  conditions.  The  fourth  group  varies 
with  the  number  of  loading  conditions.  In  this  case  the  dimension  of 
the  single  arrays  is  equal  to  the  number  of  loading  conditions.  For 
rectangular  arrays,  the  first  dimension  is  fixed,  and  the  second  dimension 
represents  the  number  of  loading  conditions.  The  fifth  group  varies  with 
the  number  of  elements.  The  sixth  group  depends  on  the  number  of  nodes. 
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The  number  of  degrees  of  freedom  determines  the  dimensions  of  the  seventh 
group.  The  number  of  degrees  of  freedom  and  the  loading  conditions 
determine  the  size  of  the  eighth  group  of  arrays.  The  first  dimension 
of  these  arrays  represents  the  degrees  of  freedom,  and  the  second 
dimension  represents  the  loading  conditions.  The  SK  matrix  in  the  last 
group  depends  on  the  number  of  degrees  of  freedom  and  the  profile  of  the 
total  structure  stiffness  matrix  which  in  turn  depends  on  the  ordering 
of  the  node  numbers  (See  the  discussion  under  subroutine  "POP"  in 
Section  5). 

The  preliminary  estimate  of  the  size  of  the  SK  array  can  be  based 
on  the  estimation  of  the  semi -bandwidth . This  would  be  an  upper  bound 
for  the  dimension  of  SK.  The  actual  dimension  of  SK  can  be  determined 
after  passing  through  the  subroutine  "POP".  This  routine  gives  a 
number  for  the  apparent  population  of  the  stiffness  matrix  from  the 
information  of  the  element  connections.  SK  must  be  dimensioned  at 
least  as  big  or  bigger  then  the  apparent  population  in  order  to  solve  the 
problem.  Usually  SK  is  the  largest  array  in  the  program,  and  its  size 
can  be  reduced  by  numbering  the  nodes  for  the  optimum  profile  of  the 
stiffness  matrix.  In  absence  of  an  adequate  procedure  for  optimization 
of  this  profile,  some  sort  of  bandwidth  optimization  is  acceptable. 

It  should  be  noted  that  the  value  of  the  variable  MAXSK  (defined  in  the 
beginning  of  the  program)  should  be  the  same  as  the  dimension  of  SK. 

When  the  dimension  of  SK  is  changed,  the  value  of  MAXSK  should  also  be 
changed. 

The  next  largest  arrays  are  FR  and  DR.  They  represent  the  applied 
force  and  the  computed  displacement  matrix  respectively.  The  dimension 


of  the  arrays  depends  on  the  number  of  degrees  of  freedom  and  the 
independent  loading  conditions.  The  first  dimension  should  be  at  least 
as  big  or  bigger  than  the  number  of  degrees  of  freedom  of  the  problem. 
Similarly  the  second  dimension  is  determined  by  the  number  of  loading 
conditions.  In  addition  the  first  dimension  should  be  the  same  as  the 
variable  NNMAX  defined  in  the  beginning  of  the  program.  Whenever  the 
dimensions  of  FR  and  DR  are  changed,  NNMAX  must  also  be  changed  accordingly. 

The  arrays  ICOL  and  IDIAG  depend  on  the  number  of  degrees  of 
freedom  of  the  problem.  Together  they  identify  the  profile  of  the 
stiffness  matrix.  For  instance,  I COL ( I ) g.ves  the  row  number  of  the  first 
non-zero  element  in  the  Ith  column  of  the  stiffness  matrix.  IDIAG(I) 

±.  L. 

gives  the  address  of  the  diagonal  element  of  the  Pr  column  of  the 
stiffness  matrix  in  the  single  array  SK. 

The  arrays  MA,  MB,  MC  and  MD  are  assigned  for  element  connections. 
NNODES  is  for  the  type  of  elements.  The  array  TH  is  for  the  sizes 
(thickness  of  plate  elements  and  cross-sectional  area  of  bars)  of  the 
elements.  The  array  MYOUNG  identifies  the  material  type  of  the  elements. 

The  remaining  arrays  are  small  and  have  minor  influence  on  the  core 
requirements. 

Frequent  Errors  Encountered  in  Using  "ANALYZE" 

1.  The  element  connections  MA,  MB,  MC  and  MD  must  be  specified 
by  starting  with  the  lowest  node  number  for  MA  and  the  next  lowest,  but 
adjacent  node  number,  for  MB.  MC  and  MD  arc  tl.cn  defined  by  continuing 
in  the  direction  established  by  MA  and  MB.  See  the  description  of  card 
sets  8,  9,  10,  and  11  in  the  input  instructions.  Section  6. 
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2.  The  boundary  degrees  of  freedom  (IBND)  must  be  in  ascending  order. 
See  the  description  of  card  set  15  in  Section  6,  Input  Instructions. 

3.  The  first  dimension  of  FR  and  DR  must  be  the  same  as  the  value 
of  the  variable  NNMAX  (defined  at  the  beginning  of  the  program). 

4.  The  value  of  MAXSK  (defined  at  the  beginning  of  the  program) 
must  be  equal  to  or  greater  than  the  value  of  the  apparent  population 
given  by  the  routine  "POP".  The  dimension  of  the  array  SK  must  be 
equal  to  the  value  given  for  MAXSK. 

5.  The  sides  of  the  shear  panels  must  be  attached  to  one  dr  more 
normal  stress  carrying  elements  such  as  posts  (bars),  membrane 
quadrilaterals  or  triangles. 
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APPENDIX  B:  LISTING  OF  THE  PROGRAM 


1 


5 


10 


IS 


70 


25 


10 


35 


40 


45 
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PPOGRAM  ANALYZE 


74/74  OPT  * 1 


FTN  4.6*446 


08/21/78  10.12.39 


PRCGRA  M ANALVZE<INPUT,  OU  TPUT  ,TAPE5=INPUT . TAPE 6* OUTPUT  I 
THE  FOLLOHING  DIM  ARE  FOR  INTERNAL  USE 

DIMENSION  AA<3,3)  ,EK<  12,121  ,Bt  12,12)  , C<12,12)  ,XI<5> ,ETA<5) ,EE <3,3) 

1 , MAA (4) ,MBH<4>  , MCC ( 4 ) , E KK ( 1 2, 12 1 , TR ANG <41 . TFR ( 10 ) , 

2 IM<10>,JN<10>,ALS<3),TITLE<8> 

rur  FOLLOWING  01 M PERTAIN  TO  THE  NUMBER  OF  MATERIALS 
DIMENSION  YOUNGM<  20) • POISON! 20 1 ,RH01 <201 

THE  FOLLOWING  DIM  ARE  THREE  TIMES  THE  NUMBER  OF  MATERIALS 
DIMENSION  ALSTRS <6  01 

THF  FOLLOWING  DIM  PERTAIN  TO  THE  NUMBER  OF  BOUNO  CONO  <NBNDRY) 
OIMFNSION  IBNOI501 

THT  FOLLOWING  DIM  PERTAIN  TO  THE  NUMBER  OF  LOADING  CONDITIONS  <L» 
DIMENSION  FOOR <12,51  ,SSX <4,5 >,SST<4,5) , SSXY <4,5 1 , SXY<5 > , KTR< 5 > . 

1 EFSTRS<5),EFFSTR<4,5>,EDR<12,5),SX<5),SY<5> ,NJLOBS<5>, 

2 ELEENG<5) ,ENGTOT<5 I , ENGSTR< 5 » , ESR<5 1 , SFTM< 5 ) 

IF  THE  number  OF  LOADING  CONDITIONS  EXCEEO  10,  THEN  CHANGE  THE 
DIMENSION  OF  EX,  EV,  EXY  in  SUBROUTINE  stress,  engg  in  subroutine 
OLSTRS  AND  TORI,  TDR2  IN  SUBROUTINE  RESTOR 

THP  following  DIM  pertain  to  the  number  of  ELEMENTS 

DIMENSION  MA<lf 0) ,MB<1 60) ,MC  < 160) ,MOI 160) ,TH< 160) .NNOOES  <1601 , 

1 M YOUNG  < 160 ) 

THE  FOLLOWING  01 M PERTAIN  TO  THE  NUMBER  OF  JOINTS 
OIMfNSION  X<90),Y<90),Z<90) 

THE  FOLLING  DIM  PERTAIN  TO  THE  NUMBER  OF  DEG  OF  FREEDOM  <NN) 
DIMENSION  IDIAG<770),ICOL<270) 

THE  FOLLOWING  OIN  PERTAIN  TO  THE  NUMBER  OF  OEG  OF  FREEDOM  <NN> 

AND  THF  NUMBER  OF  LOADING  CONDITIONS  <U 
OIMENSION  OR<270,5),FR<270,5» 

THE  FOLLOHING  DIM  PEPTAINS  TO  THE  TOTAL  STIFFNESS  MATRIX  <SKI 
OIMFNSION  SK<9110) 


THIS  PROGRAM  WAS  DEVELOPEO 

OR.  YIPPERLA  B.  VENKAYYA 

AIR  FORCE  FLIGHT  OYNAMICS  LABORATORY  < AFFOL/FBR) 
WRIGHT -PATTERSON  AIR  FORCE  BASE, DAYTON, OHIO 


INTEGFP  TYPE 

NNMAX  MUST  BF  THF  OIMENSION  OF  FR ,DR, IDI AG, ICOL 
NNMAX  * 270 

MAXSK  MUST  BE  EQUAL  OR  GREATER  THAN  THE  DIM  OF  SK 
MAXSK  = 9110 
READ  <5,21  NSTR 
1 READ  <5, 76 1 <TITLE<I1,  I = 1,81 
7F  FORMAT  <8A10> 

WRITE  <6, 77 ) <TITLE<I1<  I * 1,81 
77  F0PMAT<5X,8A10! 

KSTR*1 

®EAO <E ,2 1 MFMBS, JOINTS, NONDRY, LOADS, MM, NR, INCHES  .KIPS .NMAT.MSSTRS 
WRITE <6, 21  MEM BS, JOINTS, NBNORY.LO AOS, MM, NR, INCHES, KIPS, NMAT.MSSTRS 
READ<5,3)FEE,PMU,RHO 
IF  < RHO  , L T • . 00001  ) RHO  *0.1 
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00  77  3 3 I = l.K'AT 

ANALYZE 

59 

KX  = 3»(I-1>  ♦ 1 

ANALYZE 

60 

fri 

alotps(kx)  = f'rr, 

ANALYZE 

61 

ALS  TPS ( K X. 1 ) = 60000  . 

ANALYZE 

62 

▼ 70  7 

ALST°S(KX*2>  = 36000. 

ANALYZE 

63 

IF  C*SSTRS  .SO.  0>  GC  TO  7702 

ANALYZE 

64 

XX  = 3*N"AT 

ANALYZE 

65 

f 6 

P'AC(5,3>  (»LSTR$(I),  I = l.KXI 

ANALYZE 

66 

no  7701  i = i.kx 

ANALYZE 

67 

7781 

ALSTRS(I)  = 10C0.»ALSTSS(I> 

ANALYZE 

66 

778? 

CONTINUE 

ANALYZE 

69 

ALS  (11  = ALSTRS(l) 

ANALYZE 

70 

70 

ALSC2)  = 9 L STFS ( 2) 

ANALYZE 

71 

ALS  (3)  = ALSTKS(7) 

ANALYZE 

7? 

HPT  TO <6  , 33  3 1 ETE  .FHU.RFC. ALS(1 ) ,ALS(2) ,ALS (3» 

ANALYZE 

73 

13  3 

FORMAT  (6F1  c.  3) 

ANALYZE 

74 

TF  (NMAT  ,LE.  1>  GO  TO  7777 

ANALYZE 

75 

74 

R EA  P (r  , 3 > (YCUNGN(I) ,FCISCN)I) ,FHC1(I) , I = l.NMAT) 

ANALYZE 

76 

00  7774  I = I.NM/it 

ANALYZE 

77 

KX  = 3«(T-1 ) 4 1 

ANALYZE 

78 

“"Tins  37?)  YfUNGM (I)  .POISON)  I)  ,RH01(I)  .ALSTRS  (KX)  , 

analyze 

79 

1ALS  TPS ( K X*  1 ) ,ALSTRS(KX«2> 

ANALYZE 

60 

»0 

77A  L 

CCNTINIJP 

ANALYZE 

61 

7777 

RFAC(4,2>(NNCCFS(I> .1=1. MEMOS ) 

ANALYZE 

62 

°Ef 0(4,2) (MA( I ),I*1,MPMBS) 

ANALYZE 

63 

0 EADS,  2)  (p"(T),I  = 1,PFM8S) 

ANALYZE 

64 

RPAO(F,2)  (RC(I)  ,1=1, Mr  **9S) 

ANALYZE 

65 

•5 

pfao(c»2)  (po(ii,i=i  ,Mr*<es) 

ANALYZE 

66 

PC60(=,3) (TH(I),I=1,PFM9S) 

ANALYZE 

67 

IF  (NMAT  .LE.  1)  GO  TO  7777 

ANALYZE 

66 

CTE  A 0 (c , 2 ) (MYCLNGCI),  I = l.Mcpqj) 

ANALYZE 

69 

777  * 

00  5464  1= 1 , ME ► p$ 

ANALYZE 

90 

°0 

IF  (NMAT  a FO.  1)  MYOWC(I)  = 1 

ANALYZE 

91 

MRITF(E,33)I,6F0CES)I) ,M YCUNG ( I ) , MA ( I ) ,M9 ( I ) , NC (I ) , 70 ( I) , Th ( I ) 

ANALYZE 

9? 

6464 

CONTINUP 

ANALYZE 

93 

33 

FORMAT (7I6,4Fir ,6) 

ANALYZE 

94 

? 

F0FMAT(14Te) 

ANALYZE 

95 

04 

rgr  = crt. (10,C**6) 

ANALYZE 

96 

F = FEE 

ANALYZE 

97 

01=1.0 

ANALYZE 

96 

IF  (MM  ,LT.  3)  GO  TO  4 

ANALYZE 

99 

REA0(5,3) (X(I) ,Y(I) ,7(1) ,1=1, JOINTS) 

ANALYZE 

LOO 

iro 

GO  TO  6 

ANALYZE 

101 

4 

RE  A C (6,3)  (X(T) ,Y(I), 1 = 1, JOINTS) 

ANALYZE 

102 

00  11  1=1, JOINTS 

ANALYZE 

103 

11 

7 ( I ) =0 , 0 

ANALYZE 

104 

3 

FOPMAT  (6F10.E) 

ANALYZE 

105 

104 

6 

continue 

analyze 

106 

IF (INGRES  .FT.  1 >G0  TO  9 

ANALYZE 

107 

7000 

FORMAT  (2  OX,  3F10.4) 

ANALYZE 

106 

00  7 1=1, JOINTS 

ANALYZE 

109 

X(I)*X(I)*12,0 

ANALYZE 

110 

110 

7(T)»7)T)*12.0 

ANALYZE 

111 

7 

Y(I)=Y(I)*J2.0 

ANALYZE 

112 

9 

CONTINUE 

ANALYZE 

113 

WRITE(6, 7000)  (X(I) ,Y(T) , 7(T) , I = 1, JOINTS) 

ANALYZE 

114 

MN=MM» JOINTS 

ANALYZE 

115 

64 
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115 

NM=  NN-NBNDR  Y 

ANALYZE 

11  = 

READ (5,2  1 CI9N0CII,I=1,NBNDRY) 

ANALYZE 

117 

WPITF«6,5> 

ANALYZE 

US 

5 

FORMAT (1H1 , ///2X, 1 OH0OUNOARI £S  ///) 

ANALYZE 

119 

WRITE  (6,  10  09»  (IBNOIII  , I*l,NBNORY) 

ANALYZE 

120 

120 

00  10  1=1, NN 

ANALYZE 

121 

00  10  J*l,  LOADS 

ANALYZE 

122 

PRC  I.  J>  = 0 

analyze 

123 

10 

FR<I,JI=0 

ANALYZE 

124 

RFADC5.  21  (NJLOOS (I I ,1  = 1 , LOADS) 

ANALYZE 

125 

125 

00  21  J=  1 , LOAOS 

ANALYZE 

126 

KH=NJLOOS(J» 

ANALYZE 

127 

12 

IF  (KH-31  13,11,  14 

analyze 

126 

IS 

KX  = KH 

analyze 

123 

GO  TO  15 

ANALYZE 

130 

1 10 

14 

KX=3 

ANALYZE 

131 

15 

RFAO (5,16)  (TFRCIl.IMCI) , JMC II ,I=1,KX1 

ANALYZE 

132 

16 

FOPMAT  (3CF10. 0,21511 

ANALYZE 

133 

DO  19  1=1,  KX 

ANALYZE 

134 

KY=MM*JM(II-MM*IM(I) 

ANALYZE 

135 

115 

19 

FR(KT,  Jl  =FR(KY,JI  +TFRCI) 

ANALYZE 

136 

KH=KH-KX 

ANALYZE 

137 

IFIKHI21, 21,12 

ANALYZE 

136 

21 

CONTINUE 

ANALYZE 

139 

IFC  KIPS  «NE . 1 1 GO  TO  666 

ANALYZE 

140 

140 

00  17  1=1,  NN 

ANALYZE 

141 

DO  17  J=l,  LOAOS 

ANALYZE 

142 

17 

FR(I,JI=1000.0*FR(I,J» 

ANALYZE 

143 

666 

CONTINUE 

ANALYZE 

144 

CALL  POP  (MEMOS  , JO  I NTS  , MM,  MA  , MB  , MC  ,M0  , NNOOES  , ICOL,  IOIAG  .NONZRO.NRI 

ANALYZE 

145 

145 

I F ( NON  ZRO  ,GT.  MAXSKIGO  TO  1000 

ANALYZE 

146 

00  6 1=1  .NONZRO 

ANALYZE 

147 

a 

SKID  =0 

ANALYZE 

146 

CALL  ELSTIC  ( 1, 0, PMU, EE) 

ANALYZE 

149 

DO  120  1=1,4 

ANALYZE 

150 

150 

M AA  < I ) =1 

ANALYZE 

151 

M9BCI)*IM 

analyze 

152 

120 

MCC (11*5 

ANALYZE 

153 

MAA  (41=1 

ANALYZE 

154 

MBB  C 4 1 =4 

ANALYZE 

155 

155 

WEIGHT  = 0,0 

ANALYZE 

15= 

00  400  L = t, MEMBS 

ANALYZE 

157 

IF  ( NMAT  ,LE.  1)  GO  TO  20 

ANALYZE 

156 

KX  = M YOUNG  CL) 

ANALYZE 

159 

E = YOUNGM(KXI*10»*6 

ANALYZE 

160 

150 

PMU  = POISON(KX) 

ANALYZE 

161 

El  = E/EEF 

ANALYZE 

162 

CALL  ELS TICCEl.PMU.EE I 

ANALYZE 

163 

20 

CALL  COORD ( MAIL) »M91L) ,MCIL> ,HO (L ) ,X , Y , 2 ,AA ,XI ,ET A, AL, NNOOES (LI , 0) 

ANALYZE 

164 

I F ( NNOOES ( L 1 -3)  102,  100,124 

ANALYZE 

165 

165 

124 

CALL  OORLTL(EK,EKK,TH(L> , QUAD, MA < L ) , MB C L) , MCC U , HOC L ) • MAA , MBB , MCC, 

ANALYZE 

16= 

1X1,  ETA, NNOOES (L I ,EE,TPANG,0) 

ANALYZE 

167 

GO  TO  101 

ANALYZE 

166 

100 

CONTINUE 

ANALYZE 

169 

CALL  PLSTIFCEK  , TH (L)  . TRI ANG  , 1,2,3  , XI , ETA, EE , 0 . • 0) 

ANALYZE 

170 

170 

QUAD  = TPIANG 

ANALYZE 

171 

101 

CALL  TRNSFM(EK,AA,B,C,HM,NNOOES(L) ,12) 

ANALYZE 

172 

r 
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GO  TO  10  3 

ANALYZE 

173 

102 

CALL  ELSTIF (AA,B,C,TH(L> , MM,  AL.El) 

ANALYZE 

174 

QUAD  = AL 

ANALYZE 

175 

1/5 

10  3 

CALL  ASEMRL  (SK  ,C , MA  CL  > , MB (L 1 , MC  ( L ) . MO  <L  1 ,MM, IOI AG. NNOOES  (L  1 , 12) 

ANALYZE 

175 

.30 

FOF  MAT!/  IX  ,6E15.5/> 

ANALYZE 

177 

IF  ( NMAT  • L E.  1)  GO  TO  405 

ANALYZE 

179 

KX  = MYOUNG(L) 

analyze 

173 

RHC  = RHOl  (KX) 

ANALYZE 

160 

160 

IF(RHOKKX)  ,LF,  . 00001)  RHO=0.1 

ANALYZE 

161 

405 

HEIGHT  = HEIGHT  ♦ TH ( L ) *OUAO*RHO 

ANALYZE 

182 

400 

CONTINUE 

ANALYZE 

163 

WRITE  16, 41 0 ) HEIGHT 

ANALYZE 

164 

410 

FORMAT(tH0,10X,25HHFIGHT  OF  THE  STRUCTURE  = ,E15.5t 

ANALYZE 

165 

175 

35 

CONTINUE 

ANALYZE 

165 

C 

CALL  PRINTK(SK,IDIAG,NN> 

ANALYZE 

167 

CALL  BOUND 2 (SK ,1 BNO, NN .NBNDR Y , IDI AG, I COL) 

ANALYZE 

189 

WRITF16,  1009)  ( ICOL(I)  , 1= 1 » NM ) 

ANALYZE 

16  3 

WRITE (6,  10  09) (ID  I AG  11) ,I  = 1,NM) 

ANALYZE 

190 

190 

NONZKO=IOIAG(NM) 

ANALYZE 

191 

1009 

FOF  MAT  (IX,  10113) 

ANALYZE 

192 

nocomp=o 

ANALYZE 

193 

CALL  PEDUCE ( FR , I BND , NN, NBNORY , LOA OS,  NNMAX) 

ANALYZE 

194 

CALL  GAUSS (SK.FR, OR, ICOL , IOI AG , LOADS , NM, NNMA X .NDCOMP ) 

ANALYZE 

195 

195 

IFINOCOMP.EQ.IOI  GO  TO  2000 

analyze 

196 

CALL  RES  TOR (OR, IBND.NN, NBNORY, LOAOS, NNMAX) 

ANALYZE 

197 

CALL  RESTOR (ER, IBND.NN, NBNORY, LOAOS.  NNMAX) 

ANALYZE 

199 

DO  112  1=1, NN 

ANALYZE 

153 

00  112  J=l, LOAOS 

ANALYZE 

200 

200 

112 

OR(I,J)=OR(I,J)/EEE 

ANALYZE 

201 

DO  160  I = 1, LOAOS 

ANALYZE 

207 

ENGSTP(I)  = 0.0 

ANALYZE 

203 

00  179  J = 1 , NN 

ANALYZE 

204 

FNGSTR(I)  = ENGSTR (I ) ♦ FR(J,I)*DR(J,I) 

ANALYZE 

205 

205 

17a 

CONTINUE 

ANALYZE 

206 

ENGSTF(I)  = .5*ENGSTR(I) 

ANALYZE 

207 

180 

CONTINUE 

ANALYZE 

206 

N PA  G E = 1 

ANALYZE 

203 

LINES  = 1 

ANALYZE 

210 

210 

00  1501  1=1, LOAOS 

ANALYZE 

211 

1501 

ENGT  OT ( I ) = 0 . 

ANALYZE 

212 

DO  300  L*l, MEMOS 

ANALYZE 

213 

IF  (NMAT  . L E.  1)  GO  TO  85 

ANALYZE 

214 

KX  = MYOUNG(L) 

ANALYZE 

215 

215 

E = YOIJNGM(KX)*10**6 

ANALYZE 

216 

PMU  = POISON!  KX', 

ANALYZE 

217 

El  = E/EEE 

ANALYZE 

216 

GALL  FLS TIC (El, PMU, EE) 

ANALYZE 

219 

TYRE  = NNOOES(  L ) * 1 0 ♦ KX 

ANALYZE 

220 

220 

IF  (MSSTRS  .EO.  0)  GO  TO  86 

ANALYZE 

221 

KY  = 3MKX-1)  ♦ 1 

ANALYZE 

222 

ALS(1I  = ALSTPS(KY) 

ANALYZE 

223 

A LS ( 2 ) = ALSTRS(KY*1I 

ANALYZE 

224 

A LS  ( 3 ) = ALSTRS  ( K Y*2) 

ANALYZE 

225 

225 

GO  TO  66 

ANALYZE 

225 

65 

TVPF  = NNOOES(L) *10  ♦ 1 

ANALYZE 

227 

66 

IF< (LINES+LOAOS!  ,LT.  54  .ANO.  L .GT.  l)GO  TO  64 

ANALYZE 

229 

LINES*1 

ANALYZE 

223 
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WRITE  <6« 98) NPAGE 

ANALYZE 

230 

330 

NPAGE=NPAGE  *1 

ANALYZE 

231 

WRITE (6. 83 ) 

analyze 

232 

84 

CONTINUE 

ANALYZE 

233 

CALL  COORD (MAIL) ,MB(L>  ,MC (L) ,MO (L ) ,X , V , Z, AA , X I.ET A, AL. NNODES < LJ . 0 > 

ANALYZE 

234 

CALL  ELFORC(AA,OR,EOR. MM, MA(L) ,MB(LI . MC ( L) , MD (L) , NNOOES(L) .LOADS, 

ANALYZE 

235 

2 35 

1NNMAX) 

ANALYZE 

236 

I F (NNODES  ( L ) .LE.  31G0  TO  110 

ANALYZE 

237 

CALL  QORLTL  (EK  ,EKK,TH(LI  , QUAO.MA  ( LI  , MB1U  .MC1LI  ,MD(L1  ,MAA, M9B.MCC, 

ANALYZE 

238 

1X1,  FTA,  NNODES (L) , EE, TRANG.l ) 

ANALYZE 

233 

CALL  qlstrs(Eor,edor,xi,eta,naa,mbb,mcc,sx,sy,sxy,efstrs,£,pmu, 

ANALYZE 

240 

240 

1LOAOS.SSX.SSY,  SS X Y, EFFST R, KT R, EKK, ELE ENG, SFTM, ALS, ESR, NN ODES (L  >1 

ANALYZE 

241 

K X=KTP  (1 1 

ANALYZE 

242 

EL£ENG(1)=E LEENG1 1»»0. 5*TH(L) 

ANALYZE 

243 

ENGTOT(l)=ENGTOT(ll*ELFENG(l) 

analyze 

244 

IF  (NNODES(L)  .EQ.  51  GO  TO  220 

ANALYZE 

245 

245 

WRITE (6, 87)  L.TH(L). QUAD. TYPE, MA(LI,HB(L) , MC ( L ) , MO (L ) , SS X < KX,  1 ) , 

ANALYZE 

246 

1SSYIKX.1) ,SSXY(KX,1), (EFFSTR(I.l) , 1*1 ,4) , ELEENG ( 1 ) , SFTM< 1 I 

ANALYZE 

247 

222 

IF  (LOADS  .EQ.  1 ) GO  TO  300 

ANALYZE 

248 

DO  211  K = 2» LOADS 

ANALYZE 

249 

KX=KTR(KI 

ANALYZE 

250 

250 

ELEENG (K I =E  LEENG  (K)*0.5*TH(LI 

ANALYZE 

251 

E NGTOT (K ) = ENGTOT ( K ) *ELEENG(  K ) 

ANALYZE 

252 

IF  (NNOOES(L)  .EQ.  5)  GO  TO  225 

ANALYZE 

253 

WRITE  (6,  95)SSX(XX,K)  .SSY(KX.K)  , SSXYC  K X.K)  , CEFFSTR  ( I.K ) ,1*1,41  , 

ANALYZE 

254 

1FLEENG  (K) , SFTM (Kl 

ANALYZE 

255 

255 

GO  TO  211 

ANALYZE 

255 

225 

WPITE (6, 82 ) SSXY(KX,KI  , ( EFFSTR ( I , K I , I *1 ,4) , ELEENG ( K) , SFTM(KI 

ANALYZE 

257 

211 

CONTINUE 

ANALYZE 

250 

GO  TO  300 

analyze 

259 

220 

WRITE  (6,81)  L,TH(D  .QUAD,  TYPE,  MA(L),K6(LI  ,MC(L)  , MOIL), 

ANALYZE 

260 

260 

1SSXYJKX, 1) , (EFFSTR (1,1 1 ,1*1,4) .ELEENG (1) ,SFTM(1> 

ANALYZE 

2 61 

GO  TO  222 

ANALYZE 

262 

110 

I F ( NNOOE S < L ) .LT.  3IGO  TO  213 

ANALYZE 

263 

CALL  STRESSCE0P,XI,ETA,1,2,3,SX,SY,SXY,£FSTRS,E,PMU,ALS,ESR, 

ANALYZE 

264 

1L0A0S,ELEENG,TPIANG,3) 

ANALYZE 

265 

265 

CLFENG  ( 1 ) * ELEENG  ( 1 ) *0.  5*TH(LI 

ANALYZE 

266 

ENGTOT  <1 1 = ENGTOT  < 1 ) *ELEENG(  1 1 

ANALYZE 

267 

SFTM ( 1)  = (1.0  - ESR ( 1 ) ) / ESP ( 1 ) 

ANALYZE 

269 

WRITE  (6,  68)  L,TH(D,TPIANG,TYPE,MA(L>  ,M9  (L ) ,MC  (L ) ,SX  ( 1 ) ,SY  ( 1 ) , 

ANALYZE 

269 

1SXY(1),EFSTRS(  11  , ELEENG  (1),SFTM(  11 

ANALYZE 

270 

270 

I F ( L OA  OS  .EQ.  1 > GO  TO  300 

ANALYZE 

271 

OO  212  K*2, LOADS 

ANALYZE 

272 

SFTM(K)  * (1.0  - ESR(KII/ES9(K) 

ANALYZE 

273 

ELEENG  (Kl*  ELEENG  (K)*0.5*TH(L) 

ANALYZE 

274 

ENGTOT (K ) * ENGTOT (XI ♦ELEENG(K) 

ANALYZE 

275 

275 

212 

WRITE (6, 94) SX(KI ,SY(KI ,SX Y( K I , EFSTRS ( Kl , ELEENG ( Kl , SFTM (Kl 

ANALYZE 

275 

GO  TO  300 

ANALYZE 

277 

213 

00  215  K*l, LOADS 

ANALYZE 

2 78 

S X(KI *E* (E0R(2,KI-E0R(l,Kll/AL 

ANALYZE 

279 

ELEENG(K 1 * ( 0. 5*SX (Kl *»2/E  l*AL»TH(LI 

ANALYZE 

280 

200 

ESFIKI  * SORT ( ( SX ( K)  /ALS ( 11 1 •* 21 

ANALYZE 

281 

SFTM(K)  * (1.0  - ESR  ( Kl  I /ESR  ( K 1 

ANALYZE 

282 

IF  (SX(K)  .GE.  0.0)  GO  TO  215 

ANALYZE 

28  3 

E SF  ( K I * SORT  l (SX  (Kl/ ALS  (2)  I **2) 

ANALYZE 

784 

SFTM  (K)  = (l.o  - ESR  ( K 1 1 /ESR  (K 1 

ANALYZE 

285 

70S 

215 

FNGTOT  (Kl  *ENGTOT(K)  *ELEENG(K) 

ANALYZE 

286 

67 
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WRITF16.69)  L,TH1L),AL,TYPE,MA1L)  .MBID.SXIl)  . ELEENG 1 1 ) 

ANALYZE 

287 

IF1LOAOS  • EQ.  1 ) GO  TO  300 

ANALYZE 

288 

DO  214  K=2, LOADS 

ANALYZE 

289 

214 

WRITE  16, 931 SX1K) , ELEENG1K) .SFTM1K) 

ANALYZE 

290 

2°0 

300 

LINES*LI  Nf  S*LOAOS*l 

ANALYZE 

291 

DO  1603  KL  * 1, LOADS 

ANALYZE 

292 

1503 

WFITE(6, 1502>KL,ENGTOT(KL» .ENGSTR IKL) 

ANALYZE 

293 

1502 

FORMAT1///,  90X.39HTHE  TOTAL  ENERGY  FOR  LOAOING  CONDITION  ,I2,4H  IS 

ANALYZE 

294 

1 .E12.4,  2X,  3M1U)  ,10X,E12.4,3HIW)  > 

ANALYZE 

295 

295 

90 

LINES*l 

ANALYZE 

296 

CALL  PRNTDfi  lFR,OF,X,Y,  7.  NN  ,HH.  LOAOS,  JOINTS.  NP  AGE  , NNNAX  ) 

ANALYZE 

297 

63 

FOR  MAT(lY,4HHFM9,2X,5HTHICKt 3X.4HAREA, 2X.4HTVPE, 1X.2HMA.2X.2HMB, 

ANALYZE 

298 

12X.2HMC.7X,  2MMn,.3X,7HSIGMA-X,4X,7HSIGMA-Y,3X,8HSIGMA-XY,3X, 

ANALYZE 

299 

27HFFSTR-1 , 3X,7Mf FSTR-2, 3X, 7MEFSTR-3, 3 X, 7HEFSTR-4, 4X, 6MENERGY, 

ANALYZE 

300 

300 

36*.  2MMS) 

ANALYZE 

301 

81 

FOR  MAT  1/15  ,F7.  3,  F9. 2,51 4.22X.E  11.4,  5E  10.4,  E 10.31 

ANALYZE 

302 

82 

FORMAT(63X,E11.4.5E10.4,E10.  3) 

ANALYZE 

303 

67 

FOR  MAT  1/16,  F7.3,F9.2,5I4,3F11.A,5E10.4,E10.3> 

ANALYZE 

304 

88 

FOR  MAT  1/15,  F7.3,F9.2,4I4,4X,3E11  .4,  E 1 0 . 4, 30  X ,E1 0 .4,  E 1 0 . 31 

ANALYZE 

305 

30*5 

89 

F Of  MAT  1/15,  F7.3.F9.2, 31 4 ,8X  . E 11 . 4, 62X  , E10 . 4 ,E10.  31 

ANALYZE 

306 

9 3 

FOR  MAT  14  IX  ,E11 .4 ,62X  , E 1 0 • 4 • E 10. 3) 

ANALYZE 

307 

94 

FORMAT  I41X,  3E1 1.4,  El  0.4, 3 OX,  El  0.4,  FI  0.3) 

ANALYZE 

308 

96 

FOR  MAT  141  X,  3E11.4,5E10.4,E10.3) 

ANALYZE 

309 

96 

FORMAT 11M1 , 120X.5HPAGE  ,13/1 

ANALYZE 

310 

310 

2000 

IFIKSTR.EO.NSTR)  GO  TO  1000 

ANALYZE 

311 

KSTR*KSTR*1 

ANALYZE 

312 

GO  TO  1 

ANALYZE 

313 

1000 

CONTINUE 

ANALYZE 

314 

STOP 

ANALYZE 

315 

315 

END 

ANALYZE 

316 

.teV 
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1 

1 

SUBROUTINE  POP(NMB,JN,MM,MA,M9,MC,M0, KTYPE, 1C, IO.NZ.NR) 

POP 

2 

DIMENSION  MA(1),MB(1>  ,MC(l),MO(l>  , IC  ( 1 > , IO  (1 ) ,KT  YPE(  1 ) 

POP 

3 

IX(I, J)*I* ( J-l )+l 

POP 

6 

NZ  = 0 

POP 

5 

5 

NN*  MM*  JN 

POP 

6 

NET  = 0 

POP 

7 

DO  10  1=1,  NN 

POP 

6 

10 

IC(I)=NN 

POP 

9 

DO  50  L= 1 » MMB 

POP 

10 

in 

NNODES=2 

POP 

11 

ITFI=0 

POP 

12 

KX=IX(MM,MA(L>  ) 

POP 

13 

KY  = IX(MM,MB(L>  > 

POP 

16 

15 

IF(IC(KY)  .LT.  KX)  GO  TO  18 

POP 

15 

15 

00  19  1*1,  MM 

POP 

16 

IC(KV)=KX 

POP 

17 

19 

KY=KY*1 

POP 

18 

18 

IF(KTYPE(L)-3)  20,16,17 

POP 

19 

16 

IF( ITFI  .EO.  1 ) GO  TO  20 

POP 

20 

?0 

KY  = IX<MN,MC(L>) 

POP 

21 

ITF  1*1 

POP 

22 

NNODES=3 

POP 

23 

GO  TO  15 

POP 

26 

17 

IF1ITRI  .EO.  2 >GO  TO  20 

POP 

25 

25 

IF ( ITRI  .EO.  1 ) GO  TO  16 

POP 

25 

KY=IX(MM,MC(L>  ) 

POP 

27 

ITRI=ITRI+1 

POP 

28 

NN00ES*4 

POP 

29 

GO  TO  15 

POP 

30 

10 

16 

KY*IX(MM,MO(L>  ) 

POP 

31 

ITRI=ITRI+ 1 

POP 

32 

GO  TO  15 

POP 

33 

20 

NET* NET*  (MM**2>*( < NNODES  + ( NNODFS- 1 ) 1 / 2) 

POP 

34 

50 

CONTINUE 

POP 

35 

15 

NET=NET- (MM**2)*NR 

POP 

36 

DO  30  1*1, NN, MM 

POP 

37 

IFIICIII  .LT.  I)  GO  TO  30 

POP 

38 

K X*  I 

POP 

39 

DO  25  J*1,MM 

POP 

60 

40 

I C( KX ) *1 

POP 

41 

25 

KX*KX  + 1 

POP 

62 

10 

CONTINUE 

POP 

43 

DO  60  1*1, NN 

POP 

46 

NZ*NZ*(I-IC(I)*1) 

POP 

45 

65 

60 

IO(I >*N7 

FOP 

66 

KX*  1 NN* t NN*  1) ) /? 

POP 

67 

NET  *NET ♦ (MM*(MM*1)*JN) /2 

POP 

68 

WRITE (6,21 

POP 

49 

WRITE (6, 3)  KX.NET ,NZ 

POP 

50 

50 

WRITE  (6,  6) 

POP 

51 

WRITE  (6,51  (IC(II  ,I*t,NN» 

POP 

52 

WRITE (6, 6) 

POP 

53 

WRITE16,  51  (10(11  ,1*1, NNl 

POP 

56 

2 

FORMAT (1H1 , ////20X , 16HGR0SS  POPULATION, 4X.14HNET  POPULATION, 

POP 

55 

55 

16 X,  19MAPPARENT  POPULATION///) 

POP 

56 

3 

FOPMAK1 8X,  116,118,122//) 

POP 

57 

6 

FORMAT (//2X.16MSTARTING  ROW  NUMBERS  FOR  EACH  COLUMN///) 

POP 

56 

5 

F0RMAT(5X,  l 0112) 

POP 

59 

6 

FORMAT (//2X.61HNUMBERS  OF  01  AGONAL  ELEMENTS  IN  SINGLE  ARRAYSTIFFNE 

POP 

60 

60 

1SS  MATRIX  ///) 

POP 

61 

RETURN 

POP 

62 

ENO 

POP 

63 

BIIS  TASK  IS  BBST  QUALITY  FRACTICABU 

nu  oop  i rwMsum  w>dd,q  - — ^ 
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1 

1 

SUBROUTINE  COOPO ( K1,K2 , K3 , K4 , X , Y , Z, AA , XI ,ET A, AL ,NNO, NO) 

COORD 

2 

DIMENSION  XU)  ,Y11>  .71)1  .AA<3,3>  ,AB(3),XI(5I  , ETA (51 

COORO 

3 

XC0MP=X(K2|-X(K1> 

COORD 

4 

YC0MP=Y(K2) -Y(K1) 

COORO 

5 

5 

ZC0MP*Z(K2)*Z(K1) 

COORO 

6 

AL=SQRTt  XC0HP*»2*YC0MP**2*ZC0MP»*2> 

COORO 

Z 

AA(l,l)*XCONP/AL 

COORO 

8 

A A (1 , 21  = YCOMP/  AL 

COORO 

9 

AAI1 ,3>*ZCOMP/AL 

COORO 

10 

10 

IF(NND  .LT.  31  RET  URN 

COORO 

11 

XC0MP=X(K3>  -X(K1) 

COORD 

12 

YC0MP=Y(K31-Y(K1> 

COORO 

13 

ZC0MPx/CK3» -Z(K1) 

COORO 

14 

ALxSQPT(XCONP**2*YCOMP»*2*ZCOMP**2) 

COORO 

15 

15 

AB(1)=XC0MP/AL 

COORO 

16 

ABC  2 ) =YC OMP/AL 

COORO 

17 

AB(3)*ZC0MP/AL 

COORO 

18 

AL=SORT((  AA(1,2>«A8(3)-AA(1,3)*A8(2>  ) •• 2 ♦ ( AA ( 1 , 3) » A8 ( 1 > 

COORO 

19 

1-AA(1,1)«AB(3)  >**2*(AA (1, 1)*AB(2> -AA ( 1,2)*AB( 1) >**2) 

COORO 

20 

20 

AA(2,1>= (( AAC1.31 *»2)*A8(1)-AA(1, 1)*AA(1,3)*AB(3)-AA(1,1)» 

COORD 

21 

1AA(1«2)*AB(2)*(AA(1,2)**2)*AB(1))/AL 

COORO 

22 

AA(2,2)s  ((  AA(1,1 )**2)*AB(2)-AA(1, l)»AA(t,?)*AB(l>-AA(l,2)* 

COORO 

23 

2AA(1,3)*AB(3>*(AA(1,3>**2I*AB(2))/AL 

COORO 

24 

A A (2, 31= (( AA(1,2>»*2)*AB(3I-AA(1,2)*AA(L,3)*AB<?>-AA<1,1>* 

COORO 

25 

25 

3AA  (1 , 3>»AB  ( l»*  (AA  C 1, 1 1 **2)  « A 8 ( 311  /AL 

COORO 

26 

IFCNO  .EQ.  11RFTURN 

COORD 

27 

XI  ( 1 1 x 0.  0 

COORO 

20 

ETA ( 1 > *0 . 0 

COORO 

29 

XI(?>=(X(K2)-X(Kt>>*AA(l,n*(YtK2)-Y(Kin»AA(l,2)*(Z(K2)-Z(Kl>>*AA 

COORO 

30 

30 

1(1,3) 

COORO 

31 

ETA(2>=0.0 

COORO 

32 

XI ( 3>  = (X(K»)-X(K1> >*AA(1, 1)* (Y(K3)-Y(Kl> 1 »AA ( 1 , 2 1 ♦ ( Z (K 31 -Z (K1 ) ) • AA 

COORO 

33 

1(1,3) 

COORO 

34 

FTA(  3>  = (X(K3I-X(K1>  ) • AA  ( 2, 1 > * ( Y (K  3) -Y  (K1 1 ) *AA  (2, 2)  ♦ < Z<  K3) -Z  <K1>  > • A 

COORO 

35 

35 

1A  (2,  3) 

COORO 

36 

IFCNNO  .LE.  3) PE TURN 

COORO 

37 

XI (4)  =(X(K4  )-X(Kl ) )»AA(1 , 1) » (Y(K4) -Y (Kl> )*AA( 1,2) ♦( Z(K5) -Z(K1) )»AA 

COORO 

38 

1(1,3) 

COORO 

39 

ETA (4>  = (X(K4)-X(K1>)*AA(2,1I ♦ ( Y ( K4 > -V ( K1 ) ) »AA (2, 2 > ♦ ( Z ( K4) -Z ( K1 > > » A 

COORD 

40 

50 

1A(2,3) 

COORO 

41 

Xl(5)*(XI(2l*Xl(3)*XI(4))  /4.0 

COORO 

42 

ETA(5)=(ETA(3)*ETAC4>  )/4.0 

COORO 

43 

RETURN 

COORO 

44 

ENr 

COORO 

45 

•BIS  n«H  IS  BIST  QUALITY  PRACTICABLE 

n^KxrYnaausHHDTODoc 
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1 

1 

SUBROUTINE  ELSTIF(A,B,C,AE,HM,AL,E> 

ELSTIF 

2 

DIMENSION  A (3.3) ,8(12.12) tC 1 12. 12) 

ELSTIF 

3 

EK=AE*E/AL 

ELSTIF 

4 

00  25  1*1, MM 

ELSTIF 

5 

5 

J*I*MM 

ELSTIF 

6 

B ( 1 , I ) *EK*  A (1,1) 

ELSTIF 

7 

8(1, J)*  -8(1,1) 

ELSTIF 

S 

8(2,1)  *-8(1,1) 

ELSTIF 

3 

25 

8(2, J)  *8  (1,1) 

ELSTIF 

10 

10 

00  26  1*1, MM 

ELSTIF 

11 

DO  26  J*  1 , MM 

ELSTIF 

12 

26 

C(I,  J)*A(1  ,I)*B(1,  J) 

ELSTIF 

13 

00  36  1*1, MM 

ELSTIF 

14 

11*1 *MM 

ELSTIF 

15 

15 

DO  36  J*  1,  MM 

ELSTIF 

16 

Jl= J+MM 

ELSTIF 

IT 

C(I,J1)*-C(I,J) 

ELSTIF 

16 

C(J1,I)*-C(I,J) 

ELSTIF 

13 

36 

C(I1,J1)*C(I,J) 

ELSTIF 

20 

20 

RETURN 

ELSTIF 

21 

END 

ELSTIF 

22 

SUB  RO  j i I ME 
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1 

1 

SUBROUTINE  PLSTIF <EKK,TH,TRIANG,M A, 

NB ,MC , X, Y, EE,S HR, NONORM) 

PLSTIF 

2 

DIMENSION  EKKd2,12>  ,X(1)  ,Y(  11  ,EE(3,3)  , 

PLSTIF 

3 

1 U(6)  » A(3,3)  ,E1(3),E2<3),AX(3> 

PLSTIF 

4 

CALL  CRAMEP  (A.TRIANG,  X,Y,MA,MB,MC) 

PLSTIF 

5 

5 

00  20  I = 1*6 

PLSTIF 

6 

00  15  II  = 1*6 

PLSTIF 

7 

15 

UdI)  = 0.0 

PLSTIF 

6 

U(I>  = 1.0 

PLSTIF 

3 

El(l)  = A ( 1 , 1)  *U  ( 1 ) ♦ A<1,2)*U<3)  ♦ 

A(1,3)*U(5) 

PLSTIF 

10 

10 

E 1 ( 2)  = A(2,1)*U(2>  ♦ A(2,2)*U(4»  ♦ 

A (2, 3)*U(6) 

PLSTIF 

11 

E1C31  = All  , 1 1 *U  < 2 1 ♦ A(1,2)»U<4)  ♦ 

A ( 1 , 3)  *U(  6)  * 

A( 2, 11 *U( 11  ♦ 

PLSTIF 

12 

1 A(2,2)*U(3)  ♦ A(2,3)*UC5> 

PLSTIF 

13 

no  20  J - 1.6 

PLSTIF 

14 

DO  16  II  = 1,6 

PLSTIF 

15 

15 

16 

U(II>  = 0.0 

PLSTIF 

16 

U(J)  = 1.0 

PLSTIF 

17 

E2(  1>  = A(1,11»U(1)  ♦ A(1,2)*U(3»  ♦ 

A ( 1 , 31  *U(5 1 

PLSTIF 

10 

E 2 < 2 1 = A(2,1)*U(2>  ♦ A(2,2)*U<4)  ♦ 

A (2,3>*U<6> 

PLSTIF 

19 

E2(3)  = A(1,1>*U(2)  ♦ A (1 ,2) *U(4I  ♦ 

A ( 1, 3) *U( 6)  * 

A ( 2, 11 *U( 1)  * 

PLSTIF 

20 

20 

1 A(2,2)»U(31  ♦ A(2,3)*U<5> 

PLSTIF 

21 

EKK(I,J)  = 0.0 

PLSTIF 

22 

IF  {NONORM  .EQ.  0>  GO  TO  14 

PLSTIF 

23 

AX(1)=SHR*»2 

PLSTIF 

24 

AX(2)=2.*AX(1)  -1. 

PLSTIF 

25 

25 

AX«1I=2.*S0RT((1.-AXU»>*AX«1H 

PLSTIF 

26 

El (31  = (E  1(21 -Eld  1 l*AX  111  *E1<  3)*AX<2) 

PLSTIF 

27 

E2(3)  = (E2(2I-E2(111*AX(1I *E2(31*AXC2> 

PLSTIF 

28 

E 1 ( 1 1 = 0.0 

PLSTIF 

2) 

El  ( 21  = 0.0 

PLSTIF 

30 

30 

E2I11  * 0.0 

PLSTIF 

31 

E 2( 2 1 * 0.0 

PLSTIF 

32 

14 

00  18  K = 1,3  *. 

PLSTIF 

33 

A X ( K 1 * 0.0  - ‘ 

PLSTIF 

34 

00  17  L * 1,3 

PLSTIF 

35 

35 

17 

AXCK)  = AX(KI  * EE(K,L1*E2(L1 

PLSTIF 

36 

18 

CONTINUE 

PLSTIF 

37 

00  19  K = 1,3 

«*K,  1 * 

PLSTIF 

38 

19 

EKKd.JI  * EKK(I,JI  ♦ E1(K1*AX(KI 

PLSTIF 

39 

EKKd.JI  * EKK(I,  J)*TH*TRIANG 

PLSTIF 

40 

40 

20 

CONTINUE 

PLSTIF 

41 

00  30  I = 1,5 

% 

PLSTIF 

42 

IX  * I ♦ 1 

PLSTIF 

43 

00  30  J = IX, 6 

PLSTIF 

44 

30 

EKKIJ.II  * EKKd.JI 

PLSTIF 

45 

45 

RETURN 

PLS’IF 

46 

ENO 

PLSTIF 

47 

t 


mcnc-^ 
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1 

1 SUPPOUTINE  CRA  PER ( A, TRIANG, X , Y, HA , HB , HC) 

CRA  PER 

2 

OIPENSION  A (3,  3)  , X (1*  , Y(  1 ( 

CRA  PER 

3 

TRIANG  = 

XCHA)*CY(HB)  - T(MCI)  -T (HA ) • ( X (HB)  - X(HC>>  ♦ 

CRA  HER 

4 

1 

(XCHBI*Y(HCI  - X (HC ) * Y (MB) ) 

CRA  HER 

5 

5 A(l,l)  * 

V ( HQ ) - Y(MC) 

CRAMER 

6 

A(2,ll  * 

X ( HC)  - X (MB) 

CRA HER 

7 

A (3,1)  * 

X(MB)*YCMC>  - X(MC)*Y(MB) 

CRA HER 

8 

A ( 1 , 2)  * 

Y ( HC)  - Y(MA) 

CRAHER 

9 

AC?, 2)  * 

X ( HA)  - X (HC ) 

CRA  HER 

10 

10  A ( 3, 2)  * 

X(HC)»Y(HA)  - X ( HA) *V (HC ) 

CRAHER 

11 

A (1,3)  * 

Y ( HA ) - Y(HBC 

CRAHER 

12 

A ( 2, 3)  * 

X ( HB)  - X (HA  ) 

CRAHER 

13 

A (3 , 3 ) * 

X(HAC»Y(H8C  - X ( HB)  *Y  (HA  ) 

CRAHER 

14 

00  10  I 

= 1,3 

CRA  PER 

1$ 

15  00  10  J 

* 1.3 

CRAHER 

15 

10  A (I  , J ) = 

A ( I,  J ) /TRIANG 

CRA  PER 

17 

TRIANG  = 

(ABS(TRIANG)  ) / 2,  0 

CRAHER 

IS 

RETURN 

CRAHER 

19 

END 

CRAHER 

20 
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PAGE 

1 

t 

SUBROUTINE  OORLTL (EK .EKK, TH,  QUAD, MA, MB. MC . MO, MAA ,HBB ,MCC, XI ,ETA, 

OORLTL 

2 

INNOOES.EE, TRANG, NO) 

OORLTL 

3 

DIMENSION  EK  (12, 12),  EKK  (12,  12)  ,HAA(1>  ,M90(1I  ,HCC(1»  ,XI(5>  ,ETA(5) 

OORLTL 

4 

1,  EF  ( 3,3)  , TRANG  ( 1 ) 

qdrltl 

5 

5 

DO  125  1=1,12 

OORLTL 

6 

DO  125  J=l, 12 

OORLTL 

7 

125 

EK(I,J)=0. 

QDRLTL 

8 

NNFM=0 

OORLTL 

9 

SHF  = 1,  0 

QDRLTL 

10 

10 

I F ( NNOOES  .LE.  4)G0  TO  108 

QDRLTL 

11 

NNRM=1 

OORLTL 

12 

I F(  NNOOES  , EQ.  5)GO  TO  108 

OORLTL 

13 

IF  (NNOOES  - 7)104,105,106 

QDRLTL 

14 

104 

XCOMP*XI(3)-XI(2l 

OORLTL 

15 

15 

YC0MP=ETA(3)-ETA(2) 

OORLTL 

16 

GO  TO  107 

OORLTL 

17 

105 

XCOMP  = XI(4)-XI(3) 

OORLTL 

18 

YC0MP=ETA(4)-ETA(3> 

QDRLTL 

19 

GO  TO  107 

QDRLTL 

20 

20 

106 

yC0MP=XI(4)-XI(ll 

QDRLTL 

21 

YC0MP=ETA(4)-ETA( 1) 

QDRLTL 

22 

107 

ALL=SQRT(XCOMP**2*YCOMP»*2) 

OORLTL 

23 

SHF  = XCOMP/ALL 

QDRLTL 

24 

108 

QUAD=0. 

QDRLTL 

25 

25 

00  130  1=1,4 

OORLTL 

26 

CALL  PLSTIF(EKK,TH,TRIANG,MAA(I) ,MBB(I) ,MCC ( I ) , XI ,ET A, EE ,SHR, NNRM) 

OOR LTL 

27 

OUAO=OUAO  + TRIA  FG 

QDRLTL 

28 

TPANG(I)  =TPIANG 

QDRLTL 

29 

130 

CALL  SUM (EK ,EKK,MAA(I),MBB(I), MCC ( I ) ) 

QORLTL 

30 

30 

CALL  CONDNS(EK,FKK,MA,MB,MC.MD,NO> 

QDRLTL 

31 

RETURN 

QDRLTL 

32 

END 

QDRLTL 

33 

‘ V 
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1 1 

1 

SUBROUTINE  SUM CEK* EKK*MA ,M3*MC> 

SUM 

2 

DIMENSION  EK< 12,12) ,EKK<12*12) ,NA  (3) 

SUM 

3 

M=2 

SUM 

4 

NAC1)*2MMA-1I  + 1 

SUM 

5 

5 

NAC  2)  =2*  (MB-1)  + 1 

SUM 

6 

NA(3>*2+  CMC-1)  *1 

SUM 

7 

I H=  0 

SUM 

0 

DO  100  1=1*6 

SUM 

9 

JH-0 

SUM 

10 

10 

I F(  I «LE«  I H)  G 0 TO  30 

SUM 

11 

IHrIH*M 

SUM 

12  1 

I HHs IH/M 

SUM 

13  I 

K X=  NA  ( IHH ) 

SUM 

14 

30 

OO  90  J=  1 * 6 

SUM 

15 

15 

I F<  J .LE.  JH)GO  TO  60 

SUM 

16 

JH=  JH ♦M 

SUM 

17 

IHHsJH/M 

SUM 

16 

KY=NA ( IHH ) 

SUM 

19  | 

60 

EK<KX*KY)=EK<KX,KY)+EKK<I,J> 

SUM 

20  < 

£0 

90 

K Y=KY+1 

SUM 

21  } 

100 

KX=KX+1 

SUM 

22  * 

RETURN 

SUM 

23 

ENO 

SUM 

24  , 

1 


J 
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PAGE 

1 

1 

SU8P.OUTI NE  CONDNS IEK.EKK, MA,  H8,  HC  , MO,  NO) 

CONONS 

2 

OIHENSION  E K<  1 2, 1 2) ,EKK(12,12) 

CONONS 

3 

OO  5 1=1,8 

CONONS 

4 

00  5 J=l,8 

CONONS 

5 

5 

5 

EKK<I,  J)  = 0. 

CONONS 

6 

OET=EK(9»9)*EK(10,10>“EK<9«10)**2 

CONONS 

7 

A X=EK<9,  9) 

CONONS 

8 

EKI9,9)=EK110,10> /DET 

CONONS 

9 

E<<  1 0,  10)  = A X/OET 

CONONS 

10 

10 

EK<  9 , 1 0)  =-EKC9.10)/DET 

CONONS 

11 

EK(10,9)=EK(9,10) 

CONONS 

12 

KX=0 

CONONS 

13 

00  10  1=9,10 

CONONS 

14 

KX=KX*1 

CONONS 

15 

15 

DO  10  J=1,S 

CONONS 

16 

00  10  <=9,  10 

CONONS 

17 

10 

EKK(KX,J)=EKK(KX,J)*EK(I,K)*EK  (K,  J ) 

CONONS 

18 

I FC  NO  ,EQ.  DRETURN 

CONONS 

19 

KX=  0 

CONONS 

20 

20 

00  20  1=9,  10 

CONONS 

21 

KX=KX+1 

CONONS 

22 

00  20  J=1 , 8 

CONONS 

23 

EK<I, J)=EKK1KX,J> 

CONDNS 

24 

20 

EKMKX,J>  = 0 

CONONS 

25 

25 

00  30  1*1,8 

CONONS 

26 

00  30  J=l,8 

CONDNS 

27 

00  30  K=9, 10 

CONONS 

28 

30 

EKK  < I , J)  =EKK(  I , J)  ♦EKU.K)  *EKCK,  J) 

CONONS 

29 

00  40  1=1,8 

CONONS 

30 

30 

00  40  J=1 , 8 

CONONS 

31 

40 

EK( I , J) =EK ( X,J)*EKK(I , J) 

CONONS 

32 

I F (HC  ,Lr,  Ha)CAL*L  CHANGE (EK ,3, 5, 4, 12, 12,0) 

CONONS 

33 

IF1HO  ,LT.  H9)C4U,  OJAANGE  IEK  , 3,7,4,12,12,0) 

I F ( MO  ,LT.  MOCAUL  CT/WtatlfK .S', 7, 4.132 , 12 , 0 ) 

RETUPN  j,  4 7' 

ENO  * 

CONONS 

34 

CONONS 

35 

35 

CONONS 

36 

CONDNS 

37 
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1 

SUBROUTINE  CHANGE < EK, IX, I Y ,NNO  ,M, L, IR » 

CHANGE 

2 

OIMENSION  EKI  M,  U 

CHANGE 

3 

KX=IX 

CHANGE 

4 

KY=IY 

CHANGE 

5 

5 

M2=2*NND 

CHANGE 

9 

IF! IP  ,EO.  1) M2*L 

CHANGE 

7 

00  10  1*1,2 

CHANGE 

6 

DO  5 J=1,M2 

CHANGE 

9 

AX=EK(KX, Jl 

CHANGE 

10 

10 

EK1KX,  J»=FK(<Y,J> 

CHANGE 

11 

5 

EKtKY,  J»  =AX 

CHANGE 

12 

KX=KX*1 

CHANGE 

13 

10 

KY=KY*1 

CHANGE 

14 

IFCIR  ,EQ.  1)  RFTURN 

CHANGE 

15 

15 

KX=KX-2 

CHANGE 

16 

KY=KY-2 

CHANGE 

17 

00  20  1=1,2 

CHANGE 

18 

00  15  J=1,H2 

CHANGE 

19 

AX*EK(J,KX) 

CHANGE 

20 

20 

EK(J,KX)  *EKt  J,KV) 

CHANGE 

21 

15 

EM  J,KY)  =AX 

CHANGE 

22 

K X=K  X ♦ 1 

CHANGE 

23 

20 

KY=KY*1 

CHANGE 

24 

RETURN 

CHANGE 

25 

25 

ENC 

CHANGE 

26 
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1 

1 

SU9ROUTINE  TRNSFM CfK.AA ,B,C» MM, NND, Ml 

TRNSFM 

2 

DIMENSION  EK(12,12)  ,A*<3,3I,B(  M,  Ml  , C ( 

H,  HI 

TRNSFM 

3 

M2=2*NN0 

TRNSFM 

4 

I F ( NNO  .GT.  4)M2=8 

TRNSFM 

5 

5 

M3  = MM*NN0 

TRNSFM 

8 

IF ( NNO  .GT.  4) M3=4*MN 

TRNSFM 

7 

DO  100  1=1, M2 

TRNSFM 

5 

JA  = MM 

TRNSFM 

9 

KA=  0 

TRNSFM 

10 

10 

I A=0 

TRNSFM 

11 

DO  100  J=1,M3 

TRNSFM 

12 

B(I,J)=0.0 

TRNSFM 

13 

IF«  J-JAI  90,90,80 

TRNSFM 

14 

80 

JA= JA+HM 

TRNSFM 

15 

15 

K A=KA*2 

TRNSFM 

16 

IA=IA+MM 

TRNSFM 

17 

90 

JAA= J-IA 

TRNSFM 

IS 

00  100  K = 1 , 2 

TRNSFM 

19 

KAA  = K+KA 

TRNSFM 

20 

20 

100 

8(1 , J»=8(I , JI+EK(I,KAA)*AA(K, JAA) 

TRNSFM 

21 

DO  200  J= 1 , M3 

TRNSFM 

22 

JA=MM 

TRNSFM 

23 

KA=0 

TRNSFM 

24 

I A=0 

TRNSFM 

25 

25 

DO  200  1=1, M3 

TRNSFM 

26 

C(I,J>=0.0 

TRNSFM 

27 

IF(I-JA)  190,190,180 

TRNSFM 

28 

180 

JA= JA+MM 

TRNSFM 

29 

KA=KA+2 

TRNSFM 

30 

30 

I A=IA+MM 

TRNSFM 

31 

190 

JAf  = I-IA 

TRNSFM 

32 

00  200  K = 1 , 2 

TRNSFM 

33 

KAA  = K+KA 

TRNSFM 

34 

200 

C (I , JI=C(I , JI+AACK ,JAA)»B(KAA, J) 

TRNSFM 

35 

35 

RETURN 

TRNSFM 

36 

END  ‘ . 

TRNSFM 

37 
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1 

1 

SUBROUTINE  ASEMBL  (A. 9. HA  ,M9.HC. HO  ,HH,  ID. NNOOES. HI 

ASEMBL 

2 

DIMENSION  A (1) I0(1I.NA(4> ,NAA (31 

ASEMBL 

3 

ix(i,ji»i*(j-n*i 

ASEMBL 

4 

NNP= NNOOES 

ASEMBL 

5 

5 

I F<  NNO  .CT,  4) NN0  = 4 

ASEMBL 

6 

M2*NND*MM 

ASEMBL 

7 

NA(1I*IX(HM,MA> 

ASEMBL 

6 

NA ( 2) -IX (MM, MB) 

ASEMBL 

9 

IF( NNOOES  . GE,  3) NA( 3)»IX(MM,MC) 

ASEMBL 

10 

10 

I F( NNOOES  .GE.  4) NA(4)*IX(MN,HO) 

ASEMBL 

11 

I F ( NNOOES  .IE.  3 ) GO  TO  5 

ASEMBL 

12 

PO  4 1*1,3 

ASEMBL 

13 

KX=I/3 

ASEMBL 

14 

K Y=I / 2 

ASEMBL 

15 

15 

I F ( NA (KX  *2  I ,LT.  NA(KV*3MG0  TO  4 

ASEMBL 

16 

KH=N  A ( KX  *2  1 

ASEMBL 

17 

NA  ( KX*  2)  =NA  (KV*3) 

ASEMBL 

16 

- 

NA(KY*3)=KH 

ASEMBL 

19 

4 

CONTINUE 

ASEMBL 

20 

?0 

5 

00  10  I*  2,  NNO 

ASEMBL 

21 

10 

NAA(I-1)=NA(I)-NA(I-1(-MM 

ASEMBL 

22 

KH*MM 

ASEMBL 

23 

I AA*NA(1I 

ASEMBL 

26 

K HHs  1 

ASEMBL 

25 

?5 

TO  30  J=l,  M2 

ASEMBL 

26 

IF(J  .LE.  KH)GO  TO  15 

ASEMBL 

27 

khh=khh*i 

ASEMBL 

26 

I AA  = NA  (KHH  ) 

ASEMBL 

29 

KH=KH*MN 

ASEMBL 

30 

TO 

15 

JX*IO(IAA) -IAA*NA  ( 1)  ' 

ASEMBL 

31 

KY  = MM 

ASEMBL 

32 

00  25  1*1,  J 

ASEMBL 

33 

IF(  J .LE  .KY  .OR.  1 .LE.  KYIGO  TO  20 

ASEMBL 

34 

KX=I/MM 

ASEMBL 

35 

35 

J X=  J X*NA  A ( K X) 

ASEMBL 

36 

K V=  KY*MM 

ASEMBL 

37 

20 

A (JX)=A(  JX)  *9(1,  J)  * •’  l'K  . 

ASEMBL 

36 

25 

JX=JX*1  ’■  '■  •» 

ASEMBL 

39 

30 

IAA=IAA*1  •" 

ASEMBL 

GO 

40 

RETURN 

ASEMBL 

41 

END 

ASEMBL 

42 

§0 
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1 

SUBROUTINE  PRINTK  CSK* IDIAG.NN) 

PRINTK 

2 

OIMENSION  SKI  1 1 t IOIAG  ( 11 

PRINTK 

3 

00  SO  1=  It  NN 

PRINTK 

4 

IFII  .GT.  1»  GO  TO  65 

PRINTK 

5 

KX*1 

PRINTK 

6 

KV=1 

PRINTK 

7 

GO  TO  70 

PRINTK 

S 

65 

KX*I0IAG<I-1>*1 

PRINTK 

9 

K Y= IOIAG  < 1 1 

PRINTK 

10 

70 

MRITEI6,  3)  I 

PRINTK 

11 

SO 

WRITE (6t  2)  < SK<K>  , K = KX, KY> 

PRINTK 

12 

3 

FORMAT  <141 

PRINTK 

13 

2 

F OP  MATCIQX.IOE  12.61 

PRINTK 

16 

RETURN 

PRINTK 

15 

END 

PRINTK 

16 
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1 

SUBROUTINE  B0UN02 1 A, 18. N ,NB, 10 . IC 1 

BOUN02 

2 

DIMENSION  4(11 ,IB(1).ID(1).IC<11 

B0UN02 

3 

IH=N8 

B0UN02 

4 

NH=N 

B0UN02 

5 

00  30  JAil.NB 

B0UND2 

6 

I A=IBf  IH 1 

BOUND 2 

7 

I FI  I A .GE.  NH)  GO  TO  20 

B0UN02 

8 

KH=IA*1 

B0UN02 

9 

I F< I A ,GT.  1»  GO  TO  5 

B0UN02 

10 

K X=  1 

BOU 102 

11 

JX=1 

B0UND2 

12 

GO  TO  6 

B0UN02 

13 

5 

JX=ID(IAI-I0(IA-1» 

BOUND? 

14 

<X=IO<IA-ll  +1 

B0UN02 

15 

6 

00  10  I=KH,NH 

B0UN02 

16 

KY*  1 

BOU NO 2 

17 

IF(IC(I»  .LE.  IA1  GO  TO  7 

B0UN02 

IS 

IC<I-1)  = IC(I>-1 

BOUND  2 

19 

11  = 1 

BOUN02 

20 

KY=0 

B0UN02 

21 

GO  TO  8 

BOU  NO 2 

22 

7 

IC(I-l)  = IC(It 

80UN02 

23 

11=1-1 

BOUND? 

24 

8 

K=IC(I) 

BOU NO 2 

25 

I0«I-1»  = ID  (II-JX-KY 

BOUND? 

26 

00  10  J*K.  1 1 

BOUND? 

27 

IF«J  . EQ.  IA)  JX  = JX*1 

BOUN02 

28 

KXXsKX+JX  •-  . 

B0UN02 

29 

A ( KX ) =A<  K X X 1 

BOUND? 

30 

10 

KX=KX+1 

' * 

B0UN02 

31 

20 

NH=NH-1 

BOUND2 

32 

I H* I H— 1 

BOUND? 

33 

30 

C ONT INUE 

B0UN02 

34 

PETUFN 

END 

BOUND? 

BOUND? 

35 

36 

tlf* 
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» SUBROUTINE  REOUCE  IF, IB, N ,NB ,L, NN I 

DIMENSION  FINN.Ll  ,19111 
00  5 J»1,L 
IH=  NB 

5 NH=N 

1 I*IB  ( IM) 

IFCI-NM>  2,4,4 

2 NM1»NH-1 

00  3 K=I,NH1 
10  K1*K»1 

3 EI*,J>  *F(K 1,J» 

4 IH=IH-1 
NHsNH-1 

IF1IH.EQ.01  GO  TO  5 
15  GO  T01 

5 CONTINUE 
RETURN 
ENO 


PAGE 

1 

REOUCE 

2 

REOUCE 

3 

REDUCE 

4 

REOUCE 

5 

REOUCE 

6 

REOUCE 

7 

REOUCE 

9 

REOUCE 

3 

REOUCE 

10 

REDUCE 

11 

REOUCE 

12 

REDUCE 

13 

REOUCE 

14 

REDUCE 

IS 

REOUCE 

16 

REOUCE 

17 

REOUCE 

18 

REOUCE 

19 
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1 

1 

SUBROUTINE  GAUSS  (A.F.O.IC.I  0, l ,N ,NN»  NOCONP1 

GAUSS 

2 

DIMENSION  Alll.ICIll.  10(11.  P(NN,  U.OtNN.U 

GAUSS 

3 

I F ( NOCOMP  .EQ.  1 1 GO  TO  15 

GAUSS 

4 

00  10  1=1, N 

GAUSS 

5 

5 

11=1-1 

GAUSS 

8 

DO  9 J=  I ♦ N 

GAUSS 

7 

I F ( IC  ( Jl  .GT.  I1GO  TO  9 

GAUSS 

8 

IX=ID(  Jl  -J  + I 

GAUSS 

9 

I F(  1 1 .EQ.  OIGO  TO  8 

GAUSS 

10 

10 

00  7 K= 1, I 1 

GAUSS 

11 

IFIIC(J)  .GT.  K .OR.  IC(I)  .GT.  K1GO  TO  7 

GAUSS 

12 

KX  = IO(I)  -I  + K 

GAUSS 

13 

KY=IO(  Jl  -J4K 

GAUSS 

14 

KZ=IO(K> 

GAUSS 

15 

15 

A(IX1=A(IX)-(A(KX)*A(K71*  A(KY)I 

GAUSS 

15 

7 

CONTINUE 

GAUSS 

17 

8 

I F( I .EQ.  J 1 GO  TO  9 

GAUSS 

IS 

KZ=I0(I» 

GAUSS 

13 

I F < A (KZ(  15,6,5 

GAUSS 

20 

20 

5 

A (IX) =A( IX 1 /A  t KZ ) 

GAUSS 

21 

9 

CONTINUE 

GAUSS 

22 

10 

CONTINUE 

GAUSS 

23 

15 

DO  4 0 K«1,L 

GAUSS 

24 

DO  30  1=1, N 

GAUSS 

25 

25 

D(I,KI=F(I,K) 

GAUSS 

25 

11=1-1 

GAUSS 

27 

I F( I 1 .EQ.  01  GO  TO  30 

GAUSS 

28 

DO  20  J=1,I1 

GAUSS 

23 

IFdCdl  .GT.  J)  GO  TO  20 

GAUSS 

30 

30 

IX=I0(I)  -I*  J 

GAUSS 

31 

Q(I,K)=0(I,K)-A(IX1*0(J,  Kl 

GAUSS 

32 

20 

CONTINUE 

GAUSS 

33 

3 0 

CONTINUE 

GAUSS 

34 

40 

CONTINUE 

GAUSS 

35 

35 

00  70  1=1, N 

GAUSS 

35 

KX=ID(I1 

GAUSS 

37 

00  7 0 K=  1 , L 

GAUSS 

38 

70 

0(I,K)=0(I,K1/A(KX) 

GAUSS 

33 

DO  9 0 K=  1 , L 

GAUSS 

40 

40 

IX=N 

GAUSS 

41 

DO  90  1=2, N 

GAUSS 

42 

IX=IX-1 

GAUSS 

43 

11=1-1 

GAUSS 

44 

<x=ix 

GAUSS 

45 

45 

DO  8 0 J=1,I1 

GAUSS 

46 

KX=KX  + 1 

GAUSS 

47 

IF(IC(KX)  .GT.  IX 1 GO  TO  80 

GAUSS 

48 

<Y=IO(KX)-KX4lX 

GAUSS 

43 

D(IX,K)=D( IX,K)-A(KV1»0(KX,XI 

GAUSS 

50 

50 

80 

CONTINUE  • 

GAUSS 

51 

90 

CONTINUE  " ’ ‘ ' " ’ * * 

GAUSS 

52 

110 

RETURN  * A j , , s 

GAUSS 

53 

6 

NDCOMP=10 

GAUSS 

54 

WRITE (6, 120)1 

GAUSS 

55 

55 

120 

FORMAT (///2X, 4 6H STRUCTURE  IS  UNSTABLE,  THE  OEGREE  OF 

FREEDOM  «,I5 1 

GAUSS 

56 

RETURN 

GAUSS 

57 

END 

GAUSS 

58 

84 
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PAGE 

1 

1 

SUBROUTINE  RESTOR(  0,IB,N,NB,1,NN> 

RESTOR 

2 

DIMENSION  0<NN,U  ,I8<  1 1 , TOR  1 ( 1 0 ) , TDR2 ( 1 0 > 

RESTOR 

3 

NH=N-NB 

RESTOR 

4 

I H=  1 

RESTOR 

5 

5 1 

I = IB(IH) 

RESTOR 

S 

IF(I.GT.NH)  GO  TO  7 

RESTOR 

7 

OO  2 K=1,L 

RESTOR 

8 

TORI  < K)  =D  ( I . K) 

RESTOR 

9 

2 

0(1. 0=0. 

RESTOR 

10 

10  3 

J=I  + 1 

RESTOR 

11 

IF (J.GT.NH)  GO  TO  5 

RESTOR 

12 

00  4 K=1 ,L 

RESTOR 

13 

4 

TDR2(K)=0( J.O 

RESTOR 

14 

5 

00  6 K=1,L 

RESTOR 

15 

15 

0(J,K>  =TFP1(K) 

RESTOR 

16 

6 

TORI  (K)  = TDR2(K) 

RESTOR 

17 

IFCI.GE.NHJ  GO  TO  9 

RESTOR 

18 

1 = 1 + 1 

RESTOR 

13 

GO  TO  3 

RESTOR 

20 

20  7 

00  8 K = 1 . L 

RESTOR 

21 

8 

0 (I  . O =0  . 

RESTOR 

22 

9 

I F(  IH.GE.NB)  GO  TO  10 

RESTOR 

23 

I H=  I H ♦ 1 

RESTOR 

24 

NH=NH  + 1 

RESTOR 

25 

25 

GO  TO  1 

RESTOR 

26 

10 

CONTINUE 

RESTOR 

27 

RETURN 

RESTOR 

28 

FNO 

RESTOR 

29 

1HIS  FA8E  IS  BEST  QUALITY  PRACTICABLE 

raOtofflfVBIIIISHBDlODDC  
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PAGE 

1 

1 

SUBROUTINE  ELFORC!  AA« OR.  EORt  HM.  HA . H3. HC. NO. NNOCES. LOADS , NNI 

ELFORC 

2 

DIMENSION  AA(3, 31  .ORINN, LOAOS)  .F0RU2  .LOAOS)  .NCON  (41 

ELFORC 

3 

NCON  (11*  MM*  CHA  -11*1 

ELFORC 

4 

NCONI2)*MM»!MB  -11*1 

ELFORC 

5 

5 

I F ( NNODES  . GE.  SI  NCON!  3>  «MM*  (NC  -11*1 

ELFORC 

5 

IFINNOOES  .GE.  41  NCON  (4)  *MM*  IND  -11*1 

ELFORC 

7 

NNOzNNOOFS 

ELFORC 

8 

IFINNO  ,GT.  41 NND*4 

ELFORC 

9 

NDSP*1 

ELFORC 

10 

10 

IFINNO  .GT.  2)  N0$P*2 

ELFCRC 

11 

DO  86  K*  1.  LOAOS 

ELFORC 

12 

KH=1 

ELFORC 

13 

00  86  KK* 1 1 NNO 

ELFORC 

14 

DO  86  I* 1 » NDSP 

ELFORC 

15 

IS 

KX=NCON  IKK) 

ELFORC 

16 

EOF 1 KH.K) =0 

ELFORC 

17 

00  85  J*  1 . MM 

ELFCRC 

18 

EDF1KH,X)=E0RIKH,K>*AAII, J)*0RIKX,K> 

ELFORC 

19 

85 

K X=XX*1 

ELFORC 

20 

20  86 

KH=KM*1 

ELFORC 

21 

RETURN 

ELFORC 

22 

ENO 

ELFORC 

23 

86 
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PAGE 

1 

1 

SUBROUTI NE  STRESS ( UV, X, T , MA , MB , MC ,SX . SY, SXY, EFST , E. P. ALS .ESR. 

STRESS 

2 

1 L. ENG.TRI&NO . NNO) 

STRESS 

3 

OIMENSION  UV(12.L).X(1).Y(1>,SX(1),SY(1) ,SXY ( 1 1 , EX ( 1 0 » , EY < 1 0) , 

STRESS 

4 

IE XV  ( 1 0) , A ( 3.3)  , EFST (11 , ENG ( 1 > , ALS ( 3) ,ESR ( 1 > 

STRESS 

5 

5 

CALL  CRAMER(A,TRIANG,X,V,MA,MB,MC> 

STRESS 

6 

00  30  K=  1 , L 

STRESS 

7 

EX(K>=0. 

STRESS 

8 

EYI  K > = 0, 

STRESS 

9 

EXY ( K ) =0 • 

STRESS 

10 

10 

K X=  0 

STRESS 

11 

00  20  1=1, 3 

STRESS 

12 

IX=I  *KX 

STFESS 

13 

EX(K1=EX(K) *A( 1,I)*UV(IX,K) 

STRESS 

16 

E Y ( K ) =EY ( K 1 *A(2,I)*UV(IX*1,K) 

STRESS 

15 

15 

€XY(K)=EXY(K)*A(2fI>«UV(IX,K)*A(l,I)*UV(XX*llK) 

STRESS 

16 

20 

KX=KX*1 

STRESS 

17 

30 

CONTINUE 

STRESS 

18 

EMU=E/(1.0-P**2> 

STRESS 

19 

G=(0.5*E)/(1.0*P) 

STRESS 

20 

20 

00  40  K=  1,  L 

STRESS 

21 

SX(K)=(EX(K)+P»EY(K)  )*FHU 

STPESS 

22 

SY(K)=(»*EX(<)  *EY  (K)  >»EMU 

STRESS 

23 

40 

SXY(K)=G*EXY(K) 

STRESS 

24 

00  90  K=  1,  L 

STRESS 

25 

25 

AAX  = ALS ( 1 > 

STRESS 

26 

AAV  = ALS 1 1 ) 

STRESS 

27 

AAX Y = ALS  ( 3) 

STRESS 

28 

IF  ( SX (K ) .LT.  0.0)  AAX  = ALS(2) 

STRESS 

29 

IF  CSV ( < 1 .LT.  0.0)  AAY  = ALS(2) 

STRESS 

30 

30 

EFST(K)  = S0PT(SX(K)*»2*SY(K)*»2-SX  (K)  *S  Y ( K>  *3.«  (SX  Y(K  1 »*2)  ) 

STRESS 

31 

ENG ( K) =( SX ( K)*EX( K) *SY(K) »EV { K ) *S X Y( K > *EXY (K) )*TRIANG 

STRESS 

32 

IF(NNO  • GT • 4)ENG(K)  = (SXY(K)*EXY(K0*TRIANG 

STRESS 

33 

ESF(K)  = SORT( (SX(K)/AAX)**2  ♦ (SY (X) /AAY) **2 

STFESS 

36 

1 - ( (SX(K) *SY(K) )/(AAX»AAY) ) ♦ ( S X Y( K ) /A  AX Y ) **2) 

STRESS 

35 

35 

IF  ( NNO  .GT.  4)  ESR(K)  » ABS (SXY ( K) ) / AAXY 

STPESS 

36 

90 

CONTINUE 

STRESS 

37 

RETURN 

STRESS 

38 

END 

STRESS 

39 

IBIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
raa  oapy  nmiLU^TOocc 

4 ' ^ 
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SUBROUTINE 
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PAGE 

1 

1 

SUBROUTINE  OLSTRS ( EDR.EDOR. XI ,E TA , MAA , M39.NCC, SX , SY , SXY , EFSTRS ,E, 

QLSTRS 

2 

1PMU,L0ADS,SSX,SSY,SSXY,EFFSTR,KTR,EKK,ENG,SFTM, ALS.ESR.NND) 

QLSTRS 

3 

DIMENSION  EOR ( 12 , LOADS > , E OOR ( 1 2. L OAOS) , XI ( 1 ) , ETA ( 1 1 , MAA ( 1) , MBB< 1 > , 

QLSTRS 

4 

1MCC( 11 , S X ( l)tSY(l) ,SXV(1 > ,EFSTRS(1).SSX(4. LOADS) , SSY ( 4 , LOADS ) , 

QLSTRS 

5 

5 

2SSXY(4,LOAOS) , EFFSTR) 4 .LOADS ),KTR(1),EK<(12« 12) , ENG(l) ,ENGG(10) 

QLSTRS 

6 

3,FSP(1),SFTM(1),ALSC1) 

QLSTRS 

7 

DO  115  K=l, LOADS 

QLSTRS 

6 

SFTM(K)  = 0.0 

QLSTRS 

9 

ENG  <0=0. 

QLSTRS 

10 

10 

KX=0 

QLSTRS 

11 

DO  115  1=9,10 

QLSTRS 

12 

K X=  KX  + 1 

QLSTRS 

13 

EOF  (1,0  =0. 

QLSTRS 

14 

00  114  J= 1 , 5 

QLSTRS 

15 

15 

114 

EOR(I.K) =EOR(I,K) + EKK ( KX , J) • EOR ( J ,K> 

QLSTRS 

16 

ECF<I,0=-EDR(I,K) 

QLSTRS 

17 

115 

CONTINUE 

QLSTRS 

IS 

00  116  K=l, LOADS 

QLSTRS 

19 

EDDR(5,K)=E0R(9,K) 

QLSTRS 

20 

20 

116 

EODR  (6,K)  = FOP(  10, K) 

QLSTRS 

21 

K X=  1 

QLSTRS 

22 

KY=3 

QLSTRS 

23 

QUAD  = 0.0 

QLSTRS 

24 

DO  200  1=1,4 

QLSTRS 

25 

25 

IF(I  .LT.  4 1 GO  TO  117 

QLSTRS 

26 

KX=1 

QLSTRS 

27 

<Y=7 

QLSTRS 

28 

117 

DO  119  J= 1 , 2 

QLSTRS 

29 

DO  118  K=l, LOADS 

QLSTRS 

30 

30 

FDnR)J,K)=EOR(KX,K) 

QLSTRS 

31 

118 

F00R(J+2,K)=EDR(KY,K> 

QLSTRS 

32 

KX=KX+1 

QLSTRS 

33 

119 

KY=KY+1 

OLSTRS 

34 

CALL  STRESS(EODR,XI,ETA,HAA(II  ,M0B(II  , MCC ( 1 1 , SX ,S Y , SX Y.EFSTRS, 

QLSTRS 

35 

35 

1 E,PMU,ALS,ESR,LOAOS,ENGG,TRIANG,NNO) 

QLSTRS 

36 

QUAO  = OUAO  ♦ TP.IANG 

QLSTRS 

37 

00  201  J=l, LOADS 

QLSTRS 

38 

FNG  ( J)  =ENG  ( J)  +FN  GG  ( J) 

OLSTRS 

39 

SSX(I,J|=SX(J) 

QLSTRS 

40 

40 

SSY(I,J)=SY(J) 

OLSTRS 

41 

SSXY(I,J)=SXY(J) 

QLSTRS 

42 

EFFSTR(I,J)=EFSTRS(J1 

QLSTRS 

43 

IF(NNO  , GT . 4)EFFSTR(I,J)=ABS(SXY(J»> 

QLSTRS 

44 

SFTM(J)  = SFTM(J)  ♦ ESR< J) *T RIANG 

QLSTRS 

45 

45 

201 

CONTINUE 

QLSTRS 

46 

200 

CONTINUE 

QLSTRS 

47 

00  205  J=1 ,LOAOS 

QLSTRS 

48 

SFTM(J)  = SFTM ( J) /QUAD 

QLSTRS 

49 

SFTM(J)  = (1.0  - SFTM(JM /SFTM(  J) 

QLSTRS 

50 

50 

AMAX=0. 

QLSTRS 

51 

DO  204  1 = 1,4  ' • -.1 

QLSTRS 

52 

I F ( AMAX  ,GT.  FFFSTR(I, Jl ) GO  TO  204 

QLSTRS 

53 

AMAX  = EFFSTR(I,J> 

QLSTRS 

54 

K TM  J > *1 

QLSTRS 

55 

55 

204 

CONTINUE 

QLSTRS 

56 

205 

CONTINUE 

QLSTRS 

57 

00  210  J=l, LOADS 

QLSTRS 

58 

KX=KTR(J) 

QLSTRS 

59 

AMAX=FFFSTR (KX ,J) 

QLSTRS 

60 

50 

00  209  1*1,4 

QLSTRS 

61 

209 

EFFSTR (I, J)=EFFSTR(I, J)/AHAX 

QLSTRS 

62 

EFFSTPIKX,  J)  = AMAX 

QLSTRS 

63 

210 

CONTINUE 

QLSTRS 

64 

RETURN 

QLSTRS 

65 

55 

OLSTRS 

66 

f®ST  ® 


JBIS  15 c«®  TO  0°° 
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SUBROUTINE 
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PAGE 

1 

1 

SUBROUTINE  PRNTDR < A ,B , X, T ,Z , N, M , L , NJ , NP, NNt 

PRNTOR 

2 

OIMFNSION  A<NN,L>,8<NN,L>.X(1>,Y(1(CZ<1> 

PRNTDR 

3 

NP=  N P+ 1 

PRNTDR 

6 

LINFS=1 

PRNTDR 

5 

5 

WRITE  ( 6»  1 ) N P 

PRNTDR 

6 

WRITE  (6,  2) 

PRNTOR 

7 

DO  10  1=1, NJ 

PRNTDR 

S 

IF  <LINES*L-54>4, 3,3 

PRNTDR 

9 

3 

LINES=1 

PRNTOR 

10 

10 

WRIT  E ( 6,  1)  NP 

PRNTDR 

11 

WPITF (6,  2) 

PRNTDR 

12 

NP=NP*1 

PRNTOR 

13 

4 

KH=M*I 

PRNTDR 

14 

KHH  = KH-M*1 

PRNTDP 

15 

IS 

I F ( M ,LT.  3 )GO  TO  11 

PRNTDR 

15 

WRITE (6,  9> I,X (I) , Y ( I » ,Z f I > , ( A <J , 11 , J=KHH,KH) , ( B CJ  , 1 > , J=KHH , KH) 

PRNTOR 

17 

GO  TO  12 

PRNTDR 

16 

11 

WRITE (6,  5) I, XII) ,Y(I> , ( A(J,1),J=KHH,KH>,«  B ( J, 1 1 , J= KHH, KM) 

PRNTOR 

19 

12 

IFU  .EO.  1)  GOTO  8 

PRNTDR 

20 

20 

DO  7 K = 2 , L 

PRNTOR 

21 

IF(M  .LT.  3 ) GO  TO  13 

PRNTDR 

2? 

WRITE  (6,  6)  ( A { J,  K)  , J=KHH , KH) , ( B(J,K),  J=KHH,KH) 

PRNTOR 

23 

GO  TO  7 

PRNTDR 

24 

*3 

WRITE  (6,  15)  « A(J,K)  ,J=KHH,KH),  ( B(J,K),  J=KHH,KH) 

PRNTDR 

25 

25 

7 

CONTINUE 

PRNTDR 

26 

8 

LINES  =IINES  *L* 1 

PRNTOR 

27 

I F ( L .EO.  1)LINES  = LINES-1 

PRNTDR 

28 

10 

CONTINUE 

PRNTOR 

29 

1 

FORMAT  (1H1 , 120X,5HPAGE  ,13/) 

PRNTDR 

30 

SO 

2 

FORMAT ( IX  ,5HJOINT,8X,2H-X,8X,2H-Y,8X,2H-Z,8X,7HFORCE-X, 

PRNTOR 

31 

17X,7HFORCE-Y,7X,7HFORCE-Z.8X,7HDISPL-X,10X,7HOISPL-Y,10X, 

PRNTOR 

32 

27H0ISPL-Z//) 

PRNTOR 

33 

a 

FORMAT(/I5,F14.3,F10.3,F10.3,F12.3,F14.3,F14.3,1PE18.8, 

PRNTDR 

34 

11PF17.8,  1PE17.8) 

PRNTOR 

35 

35 

5 

FORMAT  (/I5  , F14 .3,  F 10.  3,  1 OX,  F 12,  3,  F14,  3,  14X , 1PE1 8.  8 , 1PE  17 . 8 ) 

PRNTDR 

36 

6 

FORMAT (39X ,F12.3,F14.3,F14. 3,1PE18.8, 1PE 17 . 8, 1PE17 .6 > 

PRNTDR 

37 

15 

FORMAT(39X,F12.3,F14. 3, 1 4X, 1 PE  1 6. 6 , 1PE1 7 . 6 ) 

PRNTDR 

38 

PETURN 

PRNTOR 

39 

END 

PRNTOR 

40 
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APPENDIX  C:  LISTING  OF  THE  SAMPLE  DATA 


l 


ANALYZE 

DEMO  PROBLEM-- 

INTERMEDIATE 

COMPLEXITY 

WING 

158 

OH 

3 0 

2 

.3 

111 

1 

0 

2 

1 

.5 

. 3 

. 1 

• 

60  . 

35 

• 

60. 

on 

• 

35 

• 

.ft 

. 3 

. 1 

10.5 

.3 

. 1 

3 

3 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

'» 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

4 

K 

5 

5 

5 

5 

5 

5 

5 

5 

5 

5 

5 

5 

c 

5 

5 

c 

6 

5 

5 

5 

5 

5 

5 

5 

5 

5 

c 

5 

5 

c 

c, 

5 

5 

5 

5 

5 

5 

5 

5 

ft 

6 

5 

ft 

ft 

5 

5 

5 

5 

5 

5 

5 

r. 

c 

c 

2 

2 

2 

2 

2 

2 

2 

2 

? 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

1 

7 

3 

4 

5 

6 

7 

a 

1 

2 

1 1 

t 2 

t 3 

14 

15 

16 

19 

20 

21 

22 

2 3 

24 

25 

2ft 

29 

30 

31 

32 

3 3 

34 

35 

3f 

39 

40 

41 

4? 

43 

4 4 

4ft 

6 6 

4 ) 

50 

•it 

52 

53 

54 

55 

5.6 

59 

60 

61 

62 

0 3 

64 

0 5 

06 

69  - 

20 

71 

72 

7.3 

7 4 

75 

76 

1 

3 

ft 

7 

1 

11 

13 

15 

IS 

21 

23 

25 

29 

.3  l 

33 

35 

39 

4 l 

4 3 

h5 

49 

51 

63 

55 

5 ) 

61 

6 3 

65 

*9 

71 

73 

74 

l 

19 

29 

39 

49 

5 9 

69 

5 

13 

2.3 

33 

4.3 

5.3 

6 5 

73 

o 

12 

27 

37 

47 

57 

67 

77 

1 

3 

5 

7 

9 

1 1 

13 

15 

1 7 

19 

21 

2.3 

25 

27 

20 

31 

S3 

35 

37 

5 1 

4 1 

4 1 

46 

4 7 

4 > 

C1 

5.3 

55 

57 

59 

61 

63 

65 

6 7 

65 

21 

7 3 

75 

77 

.1 

4 

c 

6 

7 

ft 

9 

10 

1 1 

12 

1 3 

16 

13 

16 

12 

1» 

21 

22 

2.3 

24 

25 

26 

27 

2 ft 

31 

.32 

! 5 

34 

15 

36 

.5  7 

3r 

41 

42 

4 3 

44 

<•5 

4b 

47 

4t 

91 

52 

53 

64 

55 

56 

57 

5ft 

61 

6? 

6 3 

64 

65 

66 

0 7 

68 

21 

72 

73 

7l. 

75 

7 6 

77 

7ft 

2 

4 

ft 

6 

2 

12 

14 

16 

20 

22 

24 

20 

3 0 

3? 

34 

3 6 

4 0 

42 

44 

46 

50 

52 

64 

5 6 

60 

62 

64 

66 

70 

72 

74 

76 

2 

20 

30 

4 0 

50 

60 

70 

6 

14 

2 4 

34 

44 

54 

Oft 

74 

10 

15 

2ft 

3P 

46 

54 

6ft 

74 

2 

4 

6 

6 

10 

12 

14 

16 

1H 

20 

2? 

24 

26 

2P 

.3  0 

32 

34 

36 

36 

4 0 

42 

44 

4 6 

4* 

50 

52 

54 

56 

53 

60 

62 

64 

6b 

60 

7 0 

22 

74 

76 

7S 

11 

12 

13 

1* 

15 

1 6 

17 

13 

21 

22 

23 

24 

25 

26 

22 

•’ft 

31 

32 

.3  3 

34 

35 

3 6 

37 

3b 

41 

42 

4 3 

44 

45 

46 

47 

4ft 

51 

52 

53 

54 

5’ 5 

56 

57 

ftp 

61 

62 

43 

64 

65 

66 

67 

6ft 

71 

72 

73 

74 

7ft 

7 6 

77 

7! 

HI 

62 

P.3 

1)4 

f 5 

ft  6 

V.  7 

84 

u 

6 

C 

10 

1 2 

14 

16 

16 

22 

24 

26 

2ft 

3? 

34 

.36 

3 0 

42 

44 

4-, 

4 8 

52 

54 

66 

56 

62 

64 

66 

6ft 

72 

74 

70 

78 

20 

30 

40 

50 

60 

7 0 

r 0 

t 4 

24 

34 

44 

54 

64 

7ft 

f . 

18 

2ft 

3ft 

4? 

5 a 

66 

7ft 

ftp. 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

(1 

0 
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APPENDIX  D : RESULTS  OF  THE  SAMPLE  PROBLEM  ANALYSIS 
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Z-DISPLACEMENTS 


LOAD  CASE  #1  LOAD  CASE  #2 


NODE  NO 

ANALYZE 

NASTRAN 

ANALYZE 

NASTRAN 

1* 

15.085 

15.170 

14.483 

14.595 

3 

16.066 

16.151 

14.752 

14.849 

5 

17.065 

17.144 

15.008 

15.086 

7 

18.133 

18.198 

15.285 

15.344 

9 

19.278 

19.328 

15.560 

15.599 

11 

15.134 

15.218 

13.960 

14.059 

13 

15.168 

15.246 

13.420 

13.500 

15 

15.213 

15.286 

12.893 

12.956 

17 

15.229 

15.297 

12.322 

12.368 

19 

10.842 

10.913 

10.710 

10.810 

21 

10.906 

10.977 

10.280 

10.368 

23 

10.914 

10.979 

9.797 

9.866 

25 

10.920 

10.982 

9.305 

9.359 

27 

10.843 

10.902 

8.738 

8.780 

29 

7.290 

7.354 

7.485 

7.578 

31 

7.353 

7.414 

7.127 

7.205 

33 

7.334 

7.386 

6.693 

6.749 

35 

7.289 

7.338 

6.244 

6.289 

37 

7.139 

7.188 

5.703 

5.738 

39 

4.434 

4.492 

4.785 

4.870 

41 

4.498 

4.547 

4.514 

4.580 

43 

4.468 

4.504 

4.158 

4.199 

45 

4.390 

4.425 

3.777 

3.810 

47 

4.183 

4.222 

3.287 

3.315 

49 

2.324 

2.378 

2.713 

2.792 

51 

2.376 

2.414 

2.512 

2.566 

53 

2.333 

2.353 

2.231 

2.256 

55 

2.242 

2.263 

1.939 

1.960 

57 

2.005 

2.034 

1.531 

1.551 

59 

.940 

.987 

1.273 

1.343 

61 

.971 

.998 

1.128 

1.169 

63 

.926 

.930 

.922 

.932 

65 

.833 

.839 

.729 

.737 

67 

.597 

.611 

.429 

.439 

69 

.150 

.175 

.327 

.368 

71 

.220 

.236 

.317 

.342 

73 

.259 

.259 

.272 

.276 

75 

.281 

.283 

.249 

.252 

77 

.187 

.195 

.125 

.130 

♦Note:  Results  are  given  for  nodes  on  the  upper  surface  only. 

The  displacement  pattern  Is  Identical  on  the  lower  surface. 

TABLE  2:  Results  From  ANALYZE  and  NASTRAN 
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